diff --git a/tools/ALARAJOYWrapper/README.md b/tools/ALARAJOYWrapper/README.md index ca2c080c..5db430c9 100644 --- a/tools/ALARAJOYWrapper/README.md +++ b/tools/ALARAJOYWrapper/README.md @@ -99,7 +99,9 @@ Either as an optional inclusion within the main ALARAJOY groupwise processing pi When the optional `-p` argument is invoked when executing `preprocess_fendl3.py`, all reactions for all nuclides written out to `cumulative_gendf_data.dsv` will be produced and saved according to the above directory structure. The highest level directory will be the same name as the `-f` TENDL data directory, with an additional "`_plots`" tag (e.g. `tendl2017/` → `tendl2017_plots/`). -To run `xs_plotting.py` as a standalone script, an input `.yaml` file must be supplied to specify the nuclides and reactions to be plotted. Additionally, the groupwise and continuous data sources to comparatively plot can be specified in this input file. These, however, are not required and will default to `cumulative_gendf_data.dsv` and `tendl2017/` respectively. The format of this input is shown below, as well as in `example_xs_plotting_input.yaml`, which can be used as a basis to supply custom plotting parameters according to their needs. Reactions are specified by their MT number, whose reference can be found at https://www.oecd-nea.org/dbdata/data/manual-endf/endf102_MT.pdf. +To run `xs_plotting.py` as a standalone script, an input `.yaml` file must be supplied to specify the nuclides and reactions to be plotted. Additionally, the groupwise and continuous data sources to comparatively plot can be specified in this input file. These, however, are not required and will default to `cumulative_gendf_data.dsv` and `tendl2017/` respectively. The format of this input is shown below, as well as in `example_xs_plotting_input.yaml`, which can be used as a basis to supply custom plotting parameters according to their needs. Reactions are specified by their MT number, whose reference can be found at https://www.oecd-nea.org/dbdata/data/manual-endf/endf102_MT.pdf. + +Relative cross-section plots can also be produced between datasets that have been produced with the same group-structure (i.e. ALARAJOY CCFE-709 cross-sections against FISPACT-II PREPRO/GROUPIE CCFE-709 cross-sections processed into an ALARAJOY DSV). To do so, when listing the DSV paths in the input YAML, put the reference data set to be divided against first. When calling `xs_plotting.py`, include the `-r` flag. ``` # Option to specify arbitrary number of groupwise data sources @@ -154,7 +156,7 @@ Using parameter `all` for any category of element, mass number, or reaction will The cross-section plotting script is run as follows: ``` -python xs_plotting.py -y /path/to/input_file.yaml +python xs_plotting.py -y /path/to/input_file.yaml -r ``` ## Application of Processed Data to ALARA Data Conversion Methods diff --git a/tools/ALARAJOYWrapper/xs_plotting.py b/tools/ALARAJOYWrapper/xs_plotting.py index b9c007ae..e327092b 100644 --- a/tools/ALARAJOYWrapper/xs_plotting.py +++ b/tools/ALARAJOYWrapper/xs_plotting.py @@ -5,6 +5,7 @@ import njoy_tools as njt import reaction_data as rxd import matplotlib.pyplot as plt +from matplotlib.patches import StepPatch from pathlib import Path def flagged_num_to_int(num): @@ -125,33 +126,47 @@ def extract_groupwise_data_from_DSV(dsv_list, KZA, MT): cross-sections were processed. emitted (str): Particle(s) emitted from the nuclear reaction corresponding to the MT number provided. + reference_group (str): Group name of the reference group against which + to divide other group structures' cross-sections. Will be an empty + string if none of the provided DSV paths in `dsv_list` contain a + '(reference)' tag in the string (i.e. example.dsv (reference)). """ groupwise_dict = {} emitted = '' + for dsv in dsv_list: with open(dsv, 'r') as f: dsv_lines = f.readlines() - group_name = dsv_lines[0].split()[-1] + group_name, processing_code, weight_function = ( + dsv_lines[0].split()[1:] + ) + _, energy_bounds = njt.load_external_group_struct(group_name) + group_name += ' (' + if processing_code != 'NJOY': + group_name += f'{processing_code}, ' + group_name += f'{weight_function} weight function)' + for line in dsv_lines[1:-1]: rxn = line.split() 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 MT == flagged_num_to_int(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' + element, A, emitted, tendl_dir, group_names, + img_ext='png', ratio_plotting=False ): """ For a given reaction's cross-section plot produced by @@ -178,6 +193,9 @@ def set_plot_save_path( img_ext (str, optional): Option to set the image filetype for the plot to be saved, limited to Matplotlib filetypes: png, ps, pdf, svg. (Defaults to 'png') + ratio_plotting (bool, optional): Option to specify the path for + plotting the ratio series between different group structures. + (Defaults to False) Returns: save_path (pathlib._local.PosixPath): Filepath for a given reaction @@ -190,27 +208,77 @@ def set_plot_save_path( filepath: CWD/tendl2017_plots/Fe/Fe56/Fe56_(n,p)_VITAMIN-J-175.png + + Likewise, the same nuclide/reaction pair plotting the cross- + section ratios between VITAMIN-J-175 and CCFE-709 groupwise + data could produce this filepath: + + CWD/tendl2017_plots/Fe/Fe56/ratio_plots/Fe56_(n,p)_ + VITAMIN-J-175_CCFE-709_ratios.png """ if isinstance(group_names, str): group_names = [group_names] + group_names = [g.replace(' ', '_').replace('/','') for g in group_names] + nuc = f'{element}{A}' nuc_dir = Path(f'{tendl_dir}_plots') / element / nuc + if ratio_plotting: + nuc_dir /= 'ratio_plots' nuc_dir.mkdir(parents=True, exist_ok=True) - return nuc_dir / f'{nuc}_(n,{emitted})_{"_".join(group_names)}.{img_ext}' + stem = str(nuc_dir / f'{nuc}_(n,{emitted})_{"_".join(group_names)}') + if ratio_plotting: + stem += '_ratios' + + return Path(stem).with_suffix(f'.{img_ext}') + +def set_plot_parameters(ax, title, ratio_plotting=False): + """ + Apply standard plotting parameters for either of the two types of plots + producable by `xs_plotting`: `plot_single_nuc_rxn_xs()` or + `plot_relative_group_xs()`. + + Arguments: + ax (matplotlib.axes._axes.Axes): Matplotlib Axes object of the plot + being constructed. + title (str): Plot title. + ratio_plotting (bool, optional): Option to produce a ratio series plot + comparing different group structures' cross-sections with + `plot_relative_group_xs()`. + (Defaults to False) + + Returns: + ax (matplotlib.axes._axes.Axes): Updated Matplotlib Axes object of the + plot being constructed. + """ + + ylabel = 'Cross-Section [b]' + if ratio_plotting: + ylabel = 'Ratio of Cross-Sections' + + if np.log10(np.ptp(ax.get_ylim())) > 1: + ax.set_yscale('log') + + ax.set_xscale('log') + ax.set_xlabel('Energy [eV]') + ax.set_ylabel(ylabel) + ax.set_title(title) + ax.grid() + ax.legend() + + return ax def plot_single_nuc_rxn_xs( ax, element, A, emitted, continuous_dict={}, groupwise_dict={} ): """ - Create and save a plot for a singular nuclide/reaction's cross-sections - vs. energy. Can be used to plot continuous TENDL data (not processed - by ALARAJOY), alongside an arbitrary number of groupwise cross- - sections according to the group structure in which they were - converted. Groupwise and continuous data can be plotted individually - if only one type is provided. + Create a plot for a singular nuclide/reaction's cross-sections vs. energy. + Can be used to plot continuous TENDL data (not processed by ALARAJOY), + alongside an arbitrary number of groupwise cross-sections according to + the group structure in which they were converted. Groupwise and + continuous data can be plotted individually if only one type is provided. Arguments: ax (matplotlib.axes._axes.Axes): Matplotlib Axes object of the plot @@ -245,9 +313,6 @@ def plot_single_nuc_rxn_xs( } (Defaults to {}) - img_ext (str, optional): Option to set the image filetype for the plot - to be saved, limited to Matplotlib filetypes: png, ps, pdf, svg. - (Defaults to 'png') Returns: ax (matplotlib.axes._axes.Axes): Updated Matplotlib Axes object of the @@ -272,15 +337,164 @@ def plot_single_nuc_rxn_xs( title += ', '.join([g for g in groupwise_dict]) + ' (Groupwise)' - ax.set_xscale('log') - ax.set_yscale('log') - ax.set_xlabel('Energy [eV]') - ax.set_ylabel('Cross-Section [b]') - ax.set_title(title) - ax.grid() - ax.legend() + return set_plot_parameters(ax, title) - return ax +def collect_all_stair_colors(ax): + """ + From a Matplotlib Axes object, determine the series colors of all + matplotlib.patches.StepPatch subplots, namely the `stair` plots for + groupwise data (as opposed to `plot` plots for continuous data). + Organize the RGBA color data into a dictionary keyed by each group + structure, with values of the tuple of 0-1 RGBA values definining the + color/transparency of each series. + + Arguments: + ax (matplotlib.axes._axes.Axes): Matplotlib Axes object of an already- + created plot with series corresponding to different group + structures' cross-section data. + + Returns: + color_dict (dict): Dictionary keyed by each group structure, with + values of the tuple of 0-1 RGBA values defining the color/ + transparency of each series. + """ + + handles, labels = ax.get_legend_handles_labels() + color_dict = dict() + for handle, label in zip(handles, labels): + if isinstance(handle, StepPatch): + color_dict[label] = handle.get_edgecolor() + + return color_dict + +def compute_groupwise_xs_ratios(groupwise_dict): + """ + Calculate the ratios of groupwise cross-sections for each group in + `groupwise_dict` that shares a group structure as a single reference + series. `ValueError` is raised if no other groupwise data matches the + group structure of the reference group structure. + + Arguments: + groupwise_dict (dict): Nested dictionary keyed at the highest level by + the name of the group structure according to which an array of + cross-sections were processed. + + Returns: + ratio_dict (dict): Nested dictionary keyed at the highest level by + the name of the group structure according to which an array of + cross-sections were processed. Similar to `groupwise_dict`, but + without a designated key for the reference group name, as all + other group names are implicitly refering to that group's cross- + section data divided by the reference group's. All data must be of + the same group structure. + reference_group (str): Group name of the reference group structure. + Will be the same as `reference_group` from Arguments if one is + supplied, otherwise, will be the group name of the first key in + `groupwise_dict`. + """ + + reference_group = next(iter(groupwise_dict)) + ref_xs, ref_energies = groupwise_dict[reference_group].values() + + ratio_dict = {} + for group_name, group_data in groupwise_dict.items(): + if group_name == reference_group: + continue + + group_xs, group_energies = group_data.values() + if np.allclose(ref_energies, group_energies, rtol=1e-5): + ratio_dict[group_name] = { + 'ratio_xs' : np.divide( + group_xs, ref_xs, + out=np.full_like(ref_xs, np.nan, dtype=float), + where=(ref_xs != 0) + )[::-1], + 'energies' : ref_energies + } + + if not ratio_dict: + raise ValueError( + 'No groupwise cross-sections provided with the same group ' \ + f'structure as the reference structure "{reference_group}".' + ) + + return ratio_dict, reference_group + +def plot_relative_group_xs( + ax, element, A, emitted, groupwise_dict, color_dict, x_limits=(None, None) +): + """ + Create a plot of the ratio series of groupwise cross-sections relative to + a reference group structure's energy-dependent cross-sections. + + Arguments: + ax (matplotlib.axes._axes.Axes): Matplotlib Axes object of the plot + being constructed. + element (str): Symbol of the element to which the nuclide being + plotted belongs. + A (str or int): Mass number for selected isonuclide. + If the target is a metastable isomer, "m" or "n" is written after + the mass number, corresponding to the first or second metastable + states. + emitted (str): Particle(s) emitted from a nuclear reaction. + continuous_dict (dict, optional): Dictionary containing an individual + nuclide's continous TENDL cross-sections and energies for a given + reaction. Formatted as: + {'xs' : continuous_xs, 'energies' : continous_energies} + + (Defaults to {}) + groupwise_dict (dict): Nested dictionary keyed at the highest level by + the name of the group structure according to which an array of + cross-sections were processed. The form of this data structure is + as follows: + { + 'group_name_1' : { + 'xs' : groupwise_xs, + 'energies' : energy_group_bounds + }, + ... + 'group_name_n' : { + 'xs' : groupwise_xs, + 'energies' : energy_group_bounds + }, + } + color_dict (dict): Dictionary keyed by each group structure, with + values of the tuple of 0-1 RGBA values defining the color/ + transparency of each series. + xlimits (tuple, optional): Option to specify the x-axis limits for + the plot. + (Defaults to (None, None)) + + Returns: + ax (matplotlib.axes._axes.Axes): Updated Matplotlib Axes object of the + plot being constructed. + """ + + if len(groupwise_dict) < 2: + raise ValueError( + 'At least two groupwise DSV files with equivalent group' \ + 'structures are required to compute a relative cross-section.' + ) + + ratio_dict, reference_group = compute_groupwise_xs_ratios(groupwise_dict) + + for group_name, group_data in ratio_dict.items(): + ax.stairs( + group_data['ratio_xs'], group_data['energies'], + baseline=None, + label=f'{group_name} / {reference_group}', + color=np.mean( + [color_dict[group_name], color_dict[reference_group]], axis=0 + ) + ) + + ax.set_xlim(x_limits) + title = ( + f'Relative Cross-Section for $^{{{A}}}${element}(n,{emitted}):\n' + f'Reference Group = {reference_group}' + ) + + return set_plot_parameters(ax, title=title, ratio_plotting=True) def check_all_tag(param): """ @@ -348,7 +562,7 @@ def find_all_mass_nums(tendl_dir, element): mass_nums.add(nuc_match.group(1)) return mass_nums - + def main(): # Only load in yaml module when executing xs_plotting.py as a script, @@ -360,6 +574,7 @@ def main(): parser = argparse.ArgumentParser() parser.add_argument('--yaml', '-y') + parser.add_argument('--ratio_plotting', '-r', action='store_true') args = parser.parse_args() with open(args.yaml, 'r') as f: @@ -414,12 +629,12 @@ def main(): tendl_dir / f'{element}{A}.tendl', flagged_num_to_int(MT) ) - groupwise_dict, emitted = extract_groupwise_data_from_DSV( - dsv_list, KZA, MT + groupwise_dict, emitted = ( + extract_groupwise_data_from_DSV(dsv_list, KZA, MT) ) if groupwise_dict: - plot_single_nuc_rxn_xs( + ax = plot_single_nuc_rxn_xs( ax, element, A, emitted, continuous_dict, groupwise_dict ) @@ -428,6 +643,20 @@ def main(): ) plt.savefig(plot_path) + if args.ratio_plotting: + ratio_fig, ratio_ax = plt.subplots(figsize=(10,6)) + color_dict = collect_all_stair_colors(ax) + plot_relative_group_xs( + ratio_ax, element, A, emitted, groupwise_dict, + color_dict, ax.get_xlim() + ) + ratio_plot_path = set_plot_save_path( + element, A, emitted, tendl_dir, + groupwise_dict.keys(), + ratio_plotting=args.ratio_plotting + ) + plt.savefig(ratio_plot_path) + print( f'Cross-section plots saved to {plot_path.parents[2]}/, ' \ 'organized by element, nuclide, reaction.'