ACEtk data generation for electron library

tools/data_library_generator/electron/README.md

@@ -0,0 +1,90 @@
+# MC/DC Electron Data Library Generator
+Converts EPRDATA14 ACE-format electron/photon/relaxation data into MC/DC's
+per-element HDF5 format for continuous-energy electron transport.
+
+## Prerequisites
+
+- Installing ACEtk from source: [link](https://github.com/njoy/ACEtk)
+- Dependencies: `pip install h5py numpy tqdm`
+- You need the EPRDATA14 library (available from [LANL Nuclear Data](https://nucleardata.lanl.gov/ace/eprdata14)).
+
+## Environment Variables
+| Variable               | Description                                        |
+|------------------------|----------------------------------------------------|
+| `MCDC_ACELIB_ELECTRON` | Path to the EPRDATA14 data file.                   |
+| `MCDC_LIB_ELECTRON`    | Path to the output directory for MC/DC HDF5 files. |
+
+## Usage
+```bash
+export MCDC_ACELIB_ELECTRON=/path/to/eprdata14/eprdata14/eprdata14
+export MCDC_LIB_ELECTRON=/path/to/mcdc/electron/library
+
+python generate.py              # Convert only missing elements
+python generate.py --rewrite    # Regenerate all files
+python generate.py --verbose    # Print detailed per-element info
+```
+
+## What it Does
+For each element (Z=1 to Z=100) in the EPRDATA14 library, the generator:
+1. Loads all elemental tables from the single concatenated EPRDATA14 file.
+2. Extracts the principal cross section energy grid and pointwise cross sections
+   (elastic, bremsstrahlung, excitation, ionization) for all reaction channels.
+3. Extracts tabulated elastic angular distributions (cosine CDFs) per incident energy (MT-528).
+4. Extracts excitation energy loss as a function of incident energy (MT-527).
+5. Extracts bremsstrahlung outgoing photon energy distributions per incident energy (MT-526).
+6. Extracts per-subshell ionization cross sections and knock-on electron
+   energy distributions (MT-534 for K-shell, MT-535+ for higher shells).
+7. Extracts atomic relaxation (fluorescence and Auger) transition data per subshell.
+8. Writes a single HDF5 file per element (e.g., `Al.h5`).
+
+## Output HDF5 Schema
+```
+<Symbol>.h5
+├── element_symbol                              (string)
+├── atomic_number                               (int)
+├── atomic_weight_ratio                         (float)
+├── electron_reactions/
+│   ├── xs_energy_grid                          (1-D array, MeV)
+│   ├── elastic_scattering/
+│   │   └── MT-528/
+│   │       ├── xs                              (1-D array, barns)
+│   │       └── angular_cosine_distribution/
+│   │           ├── energy                      (1-D array, MeV)
+│   │           ├── offset                      (1-D array, int)
+│   │           ├── cosine                      (1-D array)
+│   │           └── cdf                         (1-D array)
+│   ├── excitation/
+│   │   └── MT-527/
+│   │       ├── xs                              (1-D array, barns)
+│   │       └── energy_loss/
+│   │           ├── energy                      (1-D array, MeV)
+│   │           └── excitation_energy_loss      (1-D array, MeV)
+│   ├── bremsstrahlung/
+│   │   └── MT-526/
+│   │       ├── xs                              (1-D array, barns)
+│   │       └── energy_distribution/
+│   │           ├── energy                      (1-D array, MeV)
+│   │           ├── offset                      (1-D array, int)
+│   │           ├── energy_out                  (1-D array, MeV)
+│   │           └── cdf                         (1-D array)
+│   └── ionization/
+│       └── MT-534/  (K-shell; MT-535, MT-536, ... for higher shells)
+│           ├── xs                              (1-D array, barns)
+│           ├── binding_energy                  (float, MeV)
+│           └── energy_distribution/
+│               ├── energy                      (1-D array, MeV)
+│               ├── offset                      (1-D array, int)
+│               ├── energy_out                  (1-D array, MeV)
+│               └── cdf                         (1-D array)
+└── atomic_relaxation/
+    └── MT-534/  (K-shell; MT-535, MT-536, ... for higher shells)
+        ├── number_of_transitions               (int)
+        ├── primary_designator                  (1-D array, int)
+        ├── secondary_designator                (1-D array, int)
+        ├── energy                              (1-D array, MeV)
+        └── probability                         (1-D array)
+```
+
+## See Also
+- [Continuous Energy Theory Guide](../../docs/source/theory/cont_energy.rst)
+- [Installation Guide — CE Library Configuration](../../docs/source/install.rst)

tools/data_library_generator/electron/generate.py

@@ -0,0 +1,247 @@
+import ACEtk
+import argparse
+import h5py
+import numpy as np
+import os
+
+from tqdm import tqdm
+
+####
+
+import util
+from util import print_error, print_note
+
+parser = argparse.ArgumentParser(description="MC/DC electron data generator")
+parser.add_argument("--rewrite", dest="rewrite", action="store_true", default=False)
+parser.add_argument("--verbose", dest="verbose", action="store_true", default=False)
+args, unargs = parser.parse_known_args()
+rewrite = args.rewrite
+verbose = args.verbose
+
+# Directories
+output_dir = os.getenv("MCDC_LIB_ELECTRON")
+ace_file = os.getenv("MCDC_ACELIB_ELECTRON")
+if output_dir is None:
+    print_error("Environment variable $MCDC_LIB_ELECTRON is not set")
+if ace_file is None:
+    print_error("Environment variable $MCDC_ACELIB_ELECTRON is not set")
+# Create output directory if needed
+os.makedirs(output_dir, exist_ok=True)
+print(f"\nACE file    : {ace_file}")
+print(f"Output dir  : {output_dir}\n")
+
+# Load all tables from the concatenated EPRDATA14 file
+print("Loading EPRDATA14 tables...")
+all_tables = ACEtk.PhotoatomicTable.from_concatenated_file(ace_file)
+table_map = {t.zaid: t for t in all_tables}
+print(f"Loaded {len(table_map)} tables\n")
+
+# Select target entries
+target_entries = []
+for zaid in table_map:
+    if not zaid.endswith(".14p"):
+        continue
+
+    Z = util.decode_epr_zaid(zaid)
+    symbol = util.Z_TO_SYMBOL[Z]
+    mcdc_name = f"{symbol}.h5"
+
+    if not rewrite and os.path.exists(f"{output_dir}/{mcdc_name}"):
+        continue
+
+    target_entries.append((zaid, Z, symbol, mcdc_name))
+
+# Loop over all elements
+pbar = tqdm(
+    target_entries,
+    disable=verbose,
+    bar_format="{l_bar}{bar}| {n_fmt}/{total_fmt}{postfix}",
+)
+for zaid, Z, symbol, mcdc_name in pbar:
+
+    if not rewrite and os.path.exists(f"{output_dir}/{mcdc_name}"):
+        continue
+
+    if verbose:
+        print("\n" + "=" * 80 + "\n")
+        print(f"Create {mcdc_name} from {zaid}\n")
+    pbar.set_postfix_str(f"{mcdc_name}")
+
+    # Load ACE table
+    ace_table = table_map[zaid]
+
+    # Create MC/DC file
+    file = h5py.File(f"{output_dir}/{mcdc_name}", "w")
+
+    # ==================================================================================
+    # Basic properties
+    # ==================================================================================
+
+    header = ace_table.header
+    file.attrs["source_title"] = header.title
+    file.attrs["source_date"] = header.date
+
+    file.create_dataset("element_symbol", data=symbol)
+    file.create_dataset("atomic_number", data=Z)
+    file.create_dataset("atomic_weight_ratio", data=ace_table.atomic_weight_ratio)
+
+    # ==================================================================================
+    # Reaction groups
+    # ==================================================================================
+    # Elastic scattering : MT-528
+    # Excitation         : MT-527
+    # Bremsstrahlung     : MT-526
+    # Ionization  : MT-534 (K), MT-535 (L1), MT-536 (L2), ...
+
+    reactions = file.create_group("electron_reactions")
+
+    elastic_group = reactions.create_group("elastic_scattering")
+    excitation_group = reactions.create_group("excitation")
+    bremsstrahlung_group = reactions.create_group("bremsstrahlung")
+    ionization_group = reactions.create_group("ionization")
+
+    elastic_group.create_group("MT-528").attrs["MT"] = 528
+    excitation_group.create_group("MT-527").attrs["MT"] = 527
+    bremsstrahlung_group.create_group("MT-526").attrs["MT"] = 526
+
+    subsh_block = ace_table.electron_subshell_block
+    N_subshells = subsh_block.number_electron_subshells
+    ionization_MTs = [534 + i for i in range(N_subshells)]
+
+    for MT in ionization_MTs:
+        ionization_group.create_group(f"MT-{MT:03}").attrs["MT"] = MT
+
+    if verbose:
+        print(f"  Reaction group MTs")
+        print(f"    - Elastic scattering MT : [528]")
+        print(f"    - Excitation MT         : [527]")
+        print(f"    - Bremsstrahlung MT     : [526]")
+        print(f"    - Ionization MTs : {ionization_MTs}")
+
+    # ==================================================================================
+    # Cross sections
+    # ==================================================================================
+
+    xs0_block = ace_table.electron_cross_section_block
+
+    xs_energy = np.array(xs0_block.energies)
+    dataset = reactions.create_dataset("xs_energy_grid", data=xs_energy)
+    dataset.attrs["unit"] = "MeV"
+
+    xs = elastic_group.create_dataset("MT-528/xs", data=np.array(xs0_block.elastic))
+    xs.attrs["unit"] = "barns"
+
+    xs = excitation_group.create_dataset(
+        "MT-527/xs", data=np.array(xs0_block.excitation)
+    )
+    xs.attrs["unit"] = "barns"
+
+    xs = bremsstrahlung_group.create_dataset(
+        "MT-526/xs", data=np.array(xs0_block.bremsstrahlung)
+    )
+    xs.attrs["unit"] = "barns"
+
+    for i, MT in enumerate(ionization_MTs):
+        xs = ionization_group.create_dataset(
+            f"MT-{MT:03}/xs", data=np.array(xs0_block.electroionisation(i + 1))
+        )
+        xs.attrs["unit"] = "barns"
+
+    # ==================================================================================
+    # Elastic angular distribution
+    # ==================================================================================
+
+    angle_group = elastic_group.create_group("MT-528/angular_cosine_distribution")
+    util.load_elastic_angular_distribution(
+        ace_table.electron_elastic_angular_distribution_block, angle_group
+    )
+
+    # ==================================================================================
+    # Excitation energy loss
+    # ==================================================================================
+
+    excit_block = ace_table.electron_excitation_energy_loss_block
+    excit_group = excitation_group.create_group("MT-527/energy_loss")
+
+    dataset = excit_group.create_dataset("energy", data=np.array(excit_block.energies))
+    dataset.attrs["unit"] = "MeV"
+
+    dataset = excit_group.create_dataset(
+        "excitation_energy_loss", data=np.array(excit_block.excitation_energy_loss)
+    )
+    dataset.attrs["unit"] = "MeV"
+
+    # ==================================================================================
+    # Bremsstrahlung energy distribution
+    # ==================================================================================
+
+    brems_group = bremsstrahlung_group.create_group("MT-526/energy_distribution")
+    util.load_bremsstrahlung(
+        ace_table.bremsstrahlung_energy_distribution_block, brems_group
+    )
+
+    # ==================================================================================
+    # Ionization: binding energy and knock-on electron energy distributions
+    # ==================================================================================
+
+    for i, MT in enumerate(ionization_MTs):
+        idx = i + 1
+        MT_group = ionization_group[f"MT-{MT:03}"]
+
+        dataset = MT_group.create_dataset(
+            "binding_energy", data=subsh_block.binding_energy(idx)
+        )
+        dataset.attrs["unit"] = "MeV"
+
+        energy_dist_group = MT_group.create_group("energy_distribution")
+        util.load_electroionization_subshell(
+            ace_table.electroionisation_energy_distribution_block(idx),
+            energy_dist_group,
+        )
+
+    # ==================================================================================
+    # Atomic relaxation (subshell transition data)
+    # ==================================================================================
+
+    relax_block = ace_table.subshell_transition_data_block
+    relaxation_group = file.create_group("atomic_relaxation")
+
+    for i, MT in enumerate(ionization_MTs):
+        idx = i + 1
+        td = relax_block.transition_data(idx)
+        MT_group = relaxation_group.create_group(f"MT-{MT:03}")
+        MT_group.attrs["MT"] = MT
+
+        N_transitions = td.number_transitions
+        MT_group.create_dataset("number_of_transitions", data=N_transitions)
+
+        if N_transitions > 0:
+            primary_designators = []
+            secondary_designators = []
+            energies = []
+            probabilities = []
+            for j in range(N_transitions):
+                jdx = j + 1
+                t = td.transition(jdx)
+                primary_designators.append(td.primary_designator(jdx))
+                secondary_designators.append(td.secondary_designator(jdx))
+                energies.append(td.energy(jdx))
+                probabilities.append(td.probability(jdx))
+
+            MT_group.create_dataset(
+                "primary_designator", data=np.array(primary_designators)
+            )
+            MT_group.create_dataset(
+                "secondary_designator", data=np.array(secondary_designators)
+            )
+            dataset = MT_group.create_dataset("energy", data=np.array(energies))
+            dataset.attrs["unit"] = "MeV"
+            MT_group.create_dataset("probability", data=np.array(probabilities))
+
+    # ==================================================================================
+    # Finalize
+    # ==================================================================================
+
+    file.close()
+
+print("")