Allowing XS plotting module to handle excitation reactions

tools/ALARAJOYWrapper/README.md

@@ -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)
 
 
 

tools/ALARAJOYWrapper/preprocess_fendl3.py

@@ -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()
@@ -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():
@@ -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(
@@ -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:

tools/ALARAJOYWrapper/xs_plotting.py

@@ -6,6 +6,7 @@
 import reaction_data as rxd
 import matplotlib.pyplot as plt
 from pathlib import Path
+from openmc.data import Reaction, endf
 
 def flagged_num_to_int(num):
     """
@@ -26,14 +27,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):
     """
@@ -61,14 +63,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:
@@ -83,16 +84,37 @@ 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' : []})
-    )
-
-    return {
-        'xs'         :   xs_table['xs'],
-        'energies'   :    xs_table['E']
-    }
+    continuous_dict = {'energies' : [], 'xs' : []}
+    MT, isomeric_state = flagged_num_to_int(MT)
+    rxn = Reaction.from_endf(endf_obj, MT)
+
+    # 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[excited_state] = product
+
+        if pathways and isomeric_state < len(pathways):
+            product = pathways[list(pathways)[isomeric_state]]
+            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):
     """
@@ -141,14 +163,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'
@@ -348,7 +370,35 @@ 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):
+    """
+    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:
+        with open(dsv, 'r') as f:
+            dsv_lines = f.readlines()
+
+        for line in dsv_lines[1:-1]:
+            rxn = line.split()
+            if rxn[0] == str(pKZA):
+                MTs.add(rxn[2])
 
+    return MTs
+
 def main():
 
     # Only load in yaml module when executing xs_plotting.py as a script,
@@ -394,7 +444,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)
@@ -403,17 +453,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
                 )