Extracting KZAs directly from ENDF files, independently from filename

tools/ALARAJOYWrapper/preprocess_fendl3.py

@@ -8,8 +8,6 @@
 from pathlib import Path
 from collections import defaultdict
 
-ISOMERIC_STATES = 'mnopqrstuvwxyz'
-
 def args():
     parser = argparse.ArgumentParser()
     parser.add_argument(
@@ -48,69 +46,6 @@ def args():
     )
     return parser.parse_args()
 
-def calculate_pKZA(element, A):
-    """
-    Construct the target (parent) nuclide's KZA value. KZA values are defined
-        as ZZAAAM, where ZZ is the nuclide's atomic number, AAA is the mass
-        number, and M is the isomeric state (0 if ground state).
-
-    Arguments:
-        element (str): Chemical symbol of the target nuclide.
-        A (int or str): Mass number of the nuclide, potentially including an
-            isomeric tag, such as "m" for the first excited state, "n" for the
-            second, etc. (Note: TENDL data should only include up to a maximum
-            of the second excited state).
-
-    Returns:
-        pKZA (int): KZA value of the target nuclide.
-    """
-
-    Z = njt.elements[element]
-    A = str(A).lower()
-
-    # Metastable states classified by TENDL as m = 1, n = 2, etc.
-    # (Generally expecting only m, occasionally n, but physically,
-    # values could go higher, so isomeric_states goes up to z = 14)
-    M = 0
-    isomer_tag = next(
-        (tag for tag in ISOMERIC_STATES if tag in A), None
-    )
-    if isomer_tag:
-        M = ISOMERIC_STATES.find(isomer_tag) + 1
-
-    if M > 2:
-        warnings.warn(
-            f'Isomeric state greater than 2. Unexpected case for TENDL2017.',
-            UserWarning
-        )
-
-    A = int(A.split(isomer_tag)[0])
-    return (Z * 1000 + A) * 10 + M
-
-def interpret_KZA(kza):
-    """
-    Infer the chemical symbol and mass number from a KZA (ZZAAAM) number.
-
-    Arguments:
-        kza (int): Unique ZZAAAM for a given nuclide.
-
-    Returns:
-        element (str): Chemical symbol of the target nuclide.
-        A (str or int): Mass number for selected isotope.
-            If the target is a metastable isomer, "m" or "n" is written after 
-            the mass number, corresponding to the first or second metastable
-            states.
-    """
-
-    kza = str(kza)
-    A = kza[-4:-1]
-    Z = kza[:kza.find(A)].zfill(2)
-    element = list(njt.elements.keys())[int(Z) - 1]
-    M = int(kza[-1])
-    if M > 0:
-        A += ISOMERIC_STATES[M - 1]
-
-    return element, A
 
 def process_pendf(
     njoy_prep_input, material_id, MTs, pKZA, mt_dict, temperature, tendl_path
@@ -148,7 +83,7 @@ def process_pendf(
             is successful.
     """
 
-    element, A = interpret_KZA(pKZA)
+    element, A = tp.interpret_KZA(pKZA)
     njoy_error = ''
     njoy_input = njt.fill_input_template(
         njoy_prep_input, material_id, MTs, element, A, mt_dict, temperature
@@ -207,7 +142,7 @@ def process_gendf(
             reaction pathways for the given parent and its MTs.
     """
 
-    element, A = interpret_KZA(pKZA)
+    element, A = tp.interpret_KZA(pKZA)
     groupr_input = njt.fill_input_template(
         njoy_groupr_input, material_id, MTs, element,
         A, mt_dict, temperature, pKZA, isomer_dict
@@ -232,19 +167,19 @@ def process_gendf(
                 gendf_MTs, mt_dict, xs_line_dict, pKZA, 
                 all_rxns, eaf_nucs, isomer_dict, gendf_path
             )
-            print(f'Finished processing {element}{A}')
+            print(f'Finished processing {element}-{A}')
 
         else:
             warnings.warn(
                 f'''The requested file (MF3) is not present in the
-                ENDF file tree for {element}{A}'''
+                ENDF file tree for {element}-{A}'''
             )
             with open('mf_fail.log', 'a') as fail:
-                fail.write(f'{element}{A} \n')
+                fail.write(f'{element}-{A} \n')
 
     else:
         warnings.warn(
-            f'''Failed to convert {element}{A}.
+            f'''Failed to convert {element}-{A}.
             NJOY error message: {njoy_error}'''
         )
 
@@ -400,10 +335,7 @@ def main():
     all_rxns = defaultdict(lambda: defaultdict(dict))
 
     for file_properties in tp.search_for_files(search_dir):
-        element = file_properties['Element']
-        A = file_properties['Mass Number']
-        pKZA = calculate_pKZA(element, A)
-        endf_path = file_properties['TENDL File Path']
+        element, A, pKZA, endf_path = tuple(file_properties.values())
         TAPE20.write_bytes(endf_path.read_bytes())
 
         material_id, MTs = tp.extract_endf_specs(TAPE20)
@@ -429,7 +361,7 @@ def main():
 
         else:
             warnings.warn(
-                f'''PENDF preparation failed for {element}{A}.
+                f'''PENDF preparation failed for {element}-{A}.
                 NJOY error message: {njoy_prep_error}'''
             )
 

tools/ALARAJOYWrapper/tendl_processing.py

@@ -2,6 +2,7 @@
 import ENDFtk
 from pathlib import Path
 from reaction_data import GAS_DF
+from njoy_tools import elements
 from collections import defaultdict
 import warnings
 import numpy as np
@@ -22,27 +23,79 @@
     for key, arr in EXCITATION_DICT.items()
     for val in arr
 }
+ISOMERIC_STATES = 'mnopqrstuvwxyz'
 
-def get_isotope(stem):
+def create_endf_file_obj(path, MF):
     """
-    Extract the element name and mass number from a given filename.
+    For a TENDL (ENDF) file containing activation data for a single nuclide,
+        store a single MF's "file" data in an ENDFtk file object.
+
     Arguments:
-        stem (str): Stem of a an ENDF (TENDL) and/or PENDF file, formatted
-            as f'{element}{mass_number}.ext'
+        path (str): Filepath to an ENDF file.
+        MF (int): ENDF file number.
 
     Returns:
-        element (str): Chemical symbol of the isotope whose data is contained
-            in the file.
-        A (str): Mass number of the isotope, potentially including the letter
-            "m" at the end if the isotope is in a metastable state.
+        file (ENDFtk.tree.File): Single MF ENDFtk file object.
+        matb (int): Unique material ID extracted from the file.
     """
 
-    for i, char in enumerate(stem):
-        if char.isdigit():
-            break
+    file = None
+    tape = ENDFtk.tree.Tape.from_file(str(path))
+    matb = tape.material_numbers[0]
+    material = tape.material(matb)
+    if MF in [MF.MF for MF in material.files]:
+        file = material.file(MF)
+
+    return file, matb
+
+def calculate_KZA_from_ENDF(filepath, MF=1, MT=451):
+    """
+    Use the ZA and LISO flags from a an MF,MT section of a single nuclide ENDF
+        formatted nuclear data file to construct a unique idenfifier for the
+        nuclide in the KZA (ZZZAAAM) format.
+
+    Arguments:
+        filepath (pathlib._local.PosixPath): Path for an ENDF-formatted file.
+        MF (int, optional): Option to specify the ENDF data block ("file")
+            from which to extract identifying data. Unless specified, will
+            direct to MF = 1 ("General Information").
+            (Defaults to 1).
+        MT (int, optional): Option to specify the file section within an ENDF
+            data block (MF) from which to extract identifying data. Unless
+            specificed, will direct to MF = 451 ("Descriptive Data and
+            Directory").
+            (Defaults to 457)
+
+    Returns:
+        KZA (int): Unique ZZZAAM for a given nuclide.
+    """
+
+    endf_file_obj, _ = create_endf_file_obj(filepath, MF)
+    nuc_data = endf_file_obj.section(MT).parse()
+    return nuc_data.ZA * 10 + nuc_data.LISO
+
+def interpret_KZA(kza):
+    """
+    Infer the chemical symbol and mass number from a KZA (ZZAAAM) number.
+
+    Arguments:
+        kza (int): Unique ZZZAAAM for a given nuclide.
+
+    Returns:
+        element (str): Chemical symbol of the target nuclide.
+        A (str): Mass number for selected isotope.
+            If the target is a metastable isomer, "m" or "n" is written after 
+            the mass number, corresponding to the first or second metastable
+            states.
+    """
 
-    element = stem[:i]
-    A = stem[i:]       
+    M = kza % 10
+    za = kza // 10
+    A = str(za % 1000)
+    Z = za // 1000
+    element = list(elements.keys())[Z - 1]
+    if M > 0:
+        A += ISOMERIC_STATES[M - 1]
 
     return element, A
 
@@ -60,10 +113,11 @@ def search_for_files(dir = Path.cwd()):
     Returns:
         file_info (list of dicts): List of dictionaries containing the
             chemical symbol, mass number, and paths to the TENDL (ENDF) files
-            for a given isotope. These dictionaries are formatted as such:
+            for a given nuclide. These dictionaries are formatted as such:
                 {
-                    'Element'     : Isotope's chemical symbol,
-                    'Mass Number' : Isotope's mass number,
+                    'Element'     : Nuclide's chemical symbol,
+                    'Mass Number' : Nuclide's mass number,
+                    'pKZA'        : Nuclide's KZA value,
                     'File Paths'  : endf_path, pendf_path
                 }
     """
@@ -72,49 +126,35 @@ def search_for_files(dir = Path.cwd()):
     for suffix in ['tendl', 'endf']:
         # Iterate alphabetically for debugging to spot where process fails
         for file in sorted(dir.glob(f'*.{suffix}')):
-            element, A = get_isotope(file.stem)
+            pKZA = calculate_KZA_from_ENDF(file, 1, 451)
+            element, A = interpret_KZA(pKZA)
+
             for existing_file in file_info:
-                if (
-                    existing_file['Element'] == element
-                    and existing_file['Mass Number'] == A
-                ):
+                if existing_file['pKZA'] == pKZA:
                     warnings.warn(
-                        f'Multiple files present in {dir} for {element}-{A}'
+                        f'Multiple files present in {dir} for ' \
+                        f'{element}-{A}.'
                     )
                     break
 
             else:
+                new_filename = dir / f'{element}{A}.tendl'
+                if new_filename != file:
+                    file.rename(new_filename)
+                    warnings.warn(
+                        f'Improperly named TENDL file {file} renamed to ' \
+                        f'{new_filename} to match nuclide ID.'
+                    )
+
                 file_info.append({
-                    'Element'         :         element,
-                    'Mass Number'     :               A,
-                    'TENDL File Path' :            file
+                    'Element'         :            element,
+                    'Mass Number'     :                  A,
+                    'pKZA'            :               pKZA,
+                    'TENDL File Path' :       new_filename
                 })
 
     return file_info
 
-def create_endf_file_obj(path, MF):
-    """
-    For a TENDL (ENDF) file containing activation data for a single nuclide,
-        store a single MF's "file" data in an ENDFtk file object.
-
-    Arguments:
-        path (str): Filepath to an ENDF file.
-        MF (int): ENDF file number.
-
-    Returns:
-        file (ENDFtk.tree.File): Single MF ENDFtk file object.
-        matb (int): Unique material ID extracted from the file.
-    """
-
-    file = None
-    tape = ENDFtk.tree.Tape.from_file(str(path))
-    matb = tape.material_numbers[0]
-    material = tape.material(matb)
-    if MF in [MF.MF for MF in material.files]:
-        file = material.file(MF)
-
-    return file, matb
-
 def extract_endf_specs(path):
     """
     Extract the material ID and MT numbers from an ENDF file.