Skip to content

Allowing XS plotting module to handle excitation reactions #305

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
eitan-weinstein wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Allowing XS plotting module to handle excitation reactions #305

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
5104bc9
Adding '*' values to flagged_num_to_int() return
Jul 15, 2026
a790d3e
Fixing plotting module to handle excitation/arbitrary reactions
Jul 23, 2026
2bde4c5
Migrating multiplicity calculations to OpenMC dependency.
Jul 28, 2026
25127b6
Removing outdated endf_parserpy usage.
Jul 28, 2026
d2be63e
Undoing outdated changes to tendl_processing
Jul 28, 2026
da6280f
Streamlining null continuous_dict assignment.
Jul 28, 2026
cc005cd
Vectorizing pathways, parsing DSV with pandas
Jul 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion tools/ALARAJOYWrapper/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Dat
- Domain-specific packages
* [Endf-parserpy](https://github.com/IAEA-NDS/endf-parserpy)
* [NJOY 2016](https://github.com/njoy/NJOY2016)
* [OpenMC](https://docs.openmc.org/en/stable/quickinstall.html) (only needed if specifying a multigroup energy structure by name from the dictionary `openmc.mgxs.GROUP_STRUCTURES`)
* [OpenMC](https://docs.openmc.org/en/stable/quickinstall.html) (needed if specifying a multigroup energy structure by name from the dictionary `openmc.mgxs.GROUP_STRUCTURES` or utilizing the `xs_plotting` module)



Expand Down
7 changes: 4 additions & 3 deletions tools/ALARAJOYWrapper/preprocess_fendl3.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
from pathlib import Path
from collections import defaultdict
from subprocess import TimeoutExpired
from openmc.data import endf

def make_argparser():
parser = argparse.ArgumentParser()
Expand Down Expand Up @@ -482,6 +483,7 @@ def store_results(
dsv.write(f'{nGroups} {group_name}\n')
for parent in sorted(all_rxns):
element, A = tp.interpret_KZA(parent)
endf_obj = endf.Evaluation(tendl_dir / f'{element}{A}.tendl')
for daughter in all_rxns[parent]:
if parent != daughter:
for MT, rxn in all_rxns[parent][daughter].items():
Expand All @@ -501,8 +503,7 @@ def store_results(
)

continuous_dict = xp.extract_continuous_data(
tendl_dir / f'{element}{A}.tendl',
xp.flagged_num_to_int(MT)
endf_obj, MT
)

energies = njt.load_external_group_struct(
Expand All @@ -521,7 +522,7 @@ def store_results(
plot_path = xp.set_plot_save_path(
element, A, emitted, tendl_dir, group_name
)

plt.savefig(plot_path)

if plotting:
Expand Down
105 changes: 76 additions & 29 deletions tools/ALARAJOYWrapper/xs_plotting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,9 @@
import njoy_tools as njt
import reaction_data as rxd
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path
from openmc.data import Reaction, endf

def flagged_num_to_int(num):
"""
Expand All @@ -26,14 +28,15 @@ def flagged_num_to_int(num):
instance of '*' contained in the original value.
"""

re_match = re.match(r'^-?\d+', str(num))
num = str(num)
re_match = re.match(r'^-?\d+', num)
if not re_match:
raise ValueError(
f'Invalid flagged number {num}. Must be formatted with numeric ' \
'characters before non-numeric characters.'
)
return int(re_match.group())

return int(re_match.group()), num.count('*')

def ensure_emission_specificity(emitted, dKZA):
"""
Expand Down Expand Up @@ -61,14 +64,13 @@ def ensure_emission_specificity(emitted, dKZA):

return emitted

def extract_continuous_data(tendl_path, MT):
def extract_continuous_data(endf_obj, MT):
"""
For a given nuclide and reaction, extract its continuous-energy cross-
sections and corresponding energies from its original TENDL file.

Arguments:
tendl_path (pathlib._local.PosixPath): Path to the nuclide's original
TENDL file.
endf_obj (openmc.data.endf.Evaluation): OpenMC parsed-ENDF object.
MT (int): Reaction identifying number.

Returns:
Expand All @@ -83,16 +85,39 @@ def extract_continuous_data(tendl_path, MT):
lists.
"""

xs_table = (
tp.parse_endf_file_level_data(tendl_path)[0]
.get(MT, {})
.get('xstable', {'E' : [], 'xs' : []})
)
continuous_dict = {'energies' : [], 'xs' : []}
MT, isomeric_state = flagged_num_to_int(MT)
rxn = Reaction.from_endf(endf_obj, MT)

return {
'xs' : xs_table['xs'],
'energies' : xs_table['E']
}
# For excitation reactions, calculate specific pathway reactions by
# multiplying reaction multiplicities by MF3 cumulative cross-sections
# interpolated by the multiplicities' energy array
if isomeric_state > 0:

pathways = []
for product in rxn.products:
if product.particle not in {'neutron', 'photon', 'electron'}:
iso_flag = re.compile(r'_e(\d+)$').search(product.particle)
excited_state = int(iso_flag.group(1)) if iso_flag else 0
pathways.append((excited_state, product))

pathways.sort(key=lambda pathway: pathway[0])

@gonuke gonuke Jul 30, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wait... these excited states are not in the same order as the isomeric flags?

Then maybe the dictionary is simpler after all...

I think this is a place where a comment stating the assumptions/expectations would be warranted:

  • isomeric states are flagged in the order of their excited state value, even if the excited states don't appear in that order

@eitan-weinstein eitan-weinstein Jul 30, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rxn.products maintains the ordering of the original file, so in principle they should always be in ascending order by ENDF conventions, but I put this there purely as a safety mechanism. Mostly, I was concerned about the possibility of mismatched indexing if we're reading a deprecated MF 9/10 section that places them out of order. This is not an issue I've encountered, but given that TENDL already exhibits some weird unconventional formatting within MF 9/10 (i.e. any MT appearing in one should necessarily appear in the other, but exclusively does not), I figured guaranteeing an ascending order would probably be prudent.

@gonuke gonuke Jul 30, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that ascending order makes sense in principle, but since we are relying on this order to be consistent with the labeling of isomers, it makes me nervous that it's not very robust.

@eitan-weinstein eitan-weinstein Jul 31, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that, if this is not a reliable safeguard, should we just stick to relying on the listed order to avoid further propagating any mislabeling issues? For that, we could simplify the whole block down to something like this:

pathways = [product for product in rxn.products if product not in {'neutron', 'photon', 'electron'}]

If we're not sorting or keying based on excitation level, but rather just preserving the index, your original array suggestion may be best suited after all.

@gonuke gonuke Jul 31, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess I don't know what the correct way to do it is, and whether it's documented anywhere?

@eitan-weinstein eitan-weinstein Jul 31, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From what I can find in the ENDF-6 manual, there is no explicit ordering of LFS states listed in MF9/10, but for MF8, the subsection structure is described as follows:

The structure of each section always starts with a HEAD record and ends with a SEND record. Subsections contain data for particular final state of the reaction product (LFS). The number of subsections NS is given on the HEAD record for the section. The subsections are ordered by increasing value of LFS... (ENDF6 Manual, Section 8.2.1)

Why this is not stated explicitly in the chapters for MF9/10 is not clear, though I've only ever encountered the subsections to be increasing order of LFS. That said, this is the specific language describing subsection structuring in MF9 is as follows:

The sections are ordered by increasing MT number. Within a section for a given MT are subsections for different final states of the daughter product (LFS). (ENDF6 Manual, Section 9.1)

While there is no mention of the ordering of the subsections, below, in the "Formats" section, the LFS quantity descriptor is:

LFS Indicator to specify the level number of the nuclide (ZAP) (as defined in File 8) produced in the reaction (MT number).
LFS = 0: the final state is the ground state.
LFS = 1: the final state is the first excited state.
LFS = 2: the final state is the second excited state.
———-
———-
LFS = 98: an unspecified range of final states. (ENDF6 Manual, Section 9.2)

Again, this is not an explicit assertion of the ordering, but it seems that the ascending order in the descriptor implies a like-order of the subsections by LFS. Of course, relying on a perceived implication (even if it is consistent with the structures I've encountered in my wading through TENDL files) may not be a strong enough basis by which to make a determination of appropriateness for any ordered subsection unpacking solution.

I'll keep digging around to see if I can find other documentation in support of or countering this notion, though I'm skeptical that it would exist if not in the authoritative ENDF6 manual.


if pathways and isomeric_state < len(pathways):
product = pathways[isomeric_state][1]
energies = product.yield_.x
continuous_dict['energies'].extend(energies)
continuous_dict['xs'].extend(
product.yield_.y * rxn.xs['0K'](energies)
)

else:
mf3_xs_table = rxn.xs.get('0K')
if mf3_xs_table:
continuous_dict['energies'].extend(mf3_xs_table.x)
continuous_dict['xs'].extend(mf3_xs_table.y)

return continuous_dict

def extract_groupwise_data_from_DSV(dsv_list, KZA, MT):
"""
Expand Down Expand Up @@ -141,14 +166,14 @@ def extract_groupwise_data_from_DSV(dsv_list, KZA, MT):
dsv_pKZA, dsv_dKZA, dsv_MT, emitted = rxn[:4]
emitted = ensure_emission_specificity(emitted, dsv_dKZA)

if KZA == dsv_pKZA and MT == flagged_num_to_int(dsv_MT)[0]:
if KZA == dsv_pKZA and str(MT) == dsv_MT:
groupwise_dict[group_name] = {
'xs' : np.array(rxn[4:]).astype(float),
'energies' : energy_bounds
}
break

return groupwise_dict, ensure_emission_specificity(emitted, dsv_dKZA)
return groupwise_dict, emitted

def set_plot_save_path(
element, A, emitted, tendl_dir, group_names, img_ext='png'
Expand Down Expand Up @@ -348,7 +373,30 @@ def find_all_mass_nums(tendl_dir, element):
mass_nums.add(nuc_match.group(1))

return mass_nums


def find_all_MTs(dsv_list, pKZA):
Comment thread
eitan-weinstein marked this conversation as resolved.
"""
Given a list of preprocessed groupwise DSV files and a parent nuclide
identified by its KZA, compile all reaction identifiers (MTs) that
exist for that nuclide in any of the DSVs.

Arguments:
dsv_list (list): List of filepaths to DSV files containing
ALARAJOYWrapper-processed groupwise TENDL data.
pKZA (int): ZZZAAAM identifier of the parent nuclide.

Returns:
MTs (set): Set of all reaction types for the parent nuclide present in
any of the DSV files provided.
"""

MTs = set()
for dsv in dsv_list:
df = pd.read_csv(dsv, sep=r'\s+', skiprows=1, header=None)
MTs.update(df.loc[df[0] == pKZA, 2])

return MTs

def main():

# Only load in yaml module when executing xs_plotting.py as a script,
Expand All @@ -357,6 +405,7 @@ def main():
from yaml import safe_load

plt.rcParams.update({'figure.max_open_warning': 0})
plot_path = None

parser = argparse.ArgumentParser()
parser.add_argument('--yaml', '-y')
Expand Down Expand Up @@ -394,7 +443,7 @@ def main():

for A in mass_nums:
KZA = str((
njt.elements[element] * 1000 + flagged_num_to_int(A)
njt.elements[element] * 1000 + flagged_num_to_int(A)[0]
) * 10 + tp.ISOMERIC_STATES.find(str(A)[-1]) + 1)

MTs = adjust_dict_for_all_tag(element_dict, A)
Expand All @@ -403,17 +452,14 @@ def main():
MTs = [MTs]

if check_all_tag(MTs):
MTs = rxd.process_mt_data(rxd.load_mt_table(
njt.set_directory() / 'mt_table.csv'
)).keys()
MTs = find_all_MTs(dsv_list, KZA)

for MT in [flagged_num_to_int(MT) for MT in MTs]:
for MT in MTs:
fig, ax = plt.subplots(figsize=(10,6))

continuous_dict = extract_continuous_data(
tendl_dir / f'{element}{A}.tendl', flagged_num_to_int(MT)
endf.Evaluation(tendl_dir / f'{element}{A}.tendl'), MT
)

groupwise_dict, emitted = extract_groupwise_data_from_DSV(
dsv_list, KZA, MT
)
Expand All @@ -428,10 +474,11 @@ def main():
)
plt.savefig(plot_path)

print(
f'Cross-section plots saved to {plot_path.parents[2]}/, ' \
'organized by element, nuclide, reaction.'
)
if plot_path:
print(
f'Cross-section plots saved to {plot_path.parents[2]}/, ' \
'organized by element, nuclide, reaction.'
)


if __name__ == '__main__':
Expand Down