diff --git a/Framework/PythonInterface/plugins/algorithms/Abins.py b/Framework/PythonInterface/plugins/algorithms/Abins.py index 666f21e9785c..f8a5e3e89e70 100644 --- a/Framework/PythonInterface/plugins/algorithms/Abins.py +++ b/Framework/PythonInterface/plugins/algorithms/Abins.py @@ -7,6 +7,8 @@ from typing import Dict +import numpy as np + from mantid.api import AlgorithmFactory, FileAction, FileProperty, PythonAlgorithm, Progress from mantid.api import WorkspaceFactory, AnalysisDataService @@ -14,7 +16,7 @@ from mantid.simpleapi import ConvertUnits, GroupWorkspaces, Load from mantid.kernel import Direction, StringListValidator import abins -from abins.abinsalgorithm import AbinsAlgorithm +from abins.abinsalgorithm import AbinsAlgorithm, AtomInfo # noinspection PyPep8Naming,PyMethodMayBeStatic @@ -184,14 +186,13 @@ def PyExec(self): self.setProperty("OutputWorkspace", self._out_ws_name) prog_reporter.report("Group workspace with all required dynamical structure factors has been constructed.") - def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None): + def _fill_s_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None): """ Puts S into workspace(s). :param s_points: dynamical factor for the given atom :param workspace: workspace to be filled with S - :param protons_number: number of protons in the given type fo atom - :param nucleons_number: number of nucleons in the given type of atom + :param species: atom/isotope identity and data """ from abins.constants import FUNDAMENTALS, ONE_DIMENSIONAL_INSTRUMENTS, ONE_DIMENSIONAL_SPECTRUM @@ -201,15 +202,11 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, # only FUNDAMENTALS [data is 2d with one row] if s_points.shape[0] == FUNDAMENTALS: - self._fill_s_1d_workspace( - s_points=s_points[0], workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_1d_workspace(s_points=s_points[0], workspace=workspace, species=species) # total workspaces [data is 1d vector] elif len(s_points.shape) == ONE_DIMENSIONAL_SPECTRUM: - self._fill_s_1d_workspace( - s_points=s_points, workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_1d_workspace(s_points=s_points, workspace=workspace, species=species) # quantum order events (fundamentals or overtones + combinations for the given order) # [data is 2d table of S with a row for each quantum order] @@ -221,28 +218,19 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, wrk_name = f"{workspace}_quantum_event_{n + 1}" partial_wrk_names.append(wrk_name) - self._fill_s_1d_workspace( - s_points=s_points[n], workspace=wrk_name, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_1d_workspace(s_points=s_points[n], workspace=wrk_name, species=species) GroupWorkspaces(InputWorkspaces=partial_wrk_names, OutputWorkspace=workspace) - def _fill_s_1d_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None): + def _fill_s_1d_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None): """ Puts 1D S into workspace. - :param protons_number: number of protons in the given type of atom - :param nucleons_number: number of nucleons in the given type of atom :param s_points: dynamical factor for the given atom :param workspace: workspace to be filled with S + :param species: atom/isotope identity and data """ - if protons_number is not None: - s_points = ( - s_points - * self._scale - * self.get_cross_section( - scattering=self._scale_by_cross_section, protons_number=protons_number, nucleons_number=nucleons_number - ) - ) + if species is not None: + s_points = s_points * self._scale * self.get_cross_section(scattering=self._scale_by_cross_section, species=species) dim = 1 length = s_points.size diff --git a/Framework/PythonInterface/plugins/algorithms/Abins2D.py b/Framework/PythonInterface/plugins/algorithms/Abins2D.py index 1267a89bb73a..a2b7e791983d 100644 --- a/Framework/PythonInterface/plugins/algorithms/Abins2D.py +++ b/Framework/PythonInterface/plugins/algorithms/Abins2D.py @@ -8,6 +8,8 @@ import numbers from typing import Dict +import numpy as np + from mantid.api import AlgorithmFactory, PythonAlgorithm, Progress from mantid.api import WorkspaceFactory, AnalysisDataService from mantid.kernel import logger @@ -15,7 +17,7 @@ # noinspection PyProtectedMember from mantid.simpleapi import GroupWorkspaces import abins -from abins.abinsalgorithm import AbinsAlgorithm +from abins.abinsalgorithm import AbinsAlgorithm, AtomInfo # noinspection PyPep8Naming,PyMethodMayBeStatic @@ -184,14 +186,13 @@ def PyExec(self): self.setProperty("OutputWorkspace", self._out_ws_name) prog_reporter.report("Group workspace with all required dynamical structure factors has been constructed.") - def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None): + def _fill_s_workspace(self, *, s_points=None, workspace: str, species: AtomInfo | None = None): """ Puts S into workspace(s). :param s_points: dynamical factor for the given atom :param workspace: workspace to be filled with S - :param protons_number: number of protons in the given type fo atom - :param nucleons_number: number of nucleons in the given type of atom + :param species: Atom/isotope data """ from abins.constants import FUNDAMENTALS, TWO_DIMENSIONAL_INSTRUMENTS @@ -200,15 +201,11 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, # only FUNDAMENTALS [data is 3d with length 1 in axis 0] if s_points.shape[0] == FUNDAMENTALS: - self._fill_s_2d_workspace( - s_points=s_points[0], workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_2d_workspace(s_points=s_points[0], workspace=workspace, species=species) # total workspaces [data is 2d array of S] elif s_points.shape[0] == abins.parameters.instruments[self._instrument.get_name()]["q_size"]: - self._fill_s_2d_workspace( - s_points=s_points, workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_2d_workspace(s_points=s_points, workspace=workspace, species=species) # Multiple quantum order events [data is 3d table of S using axis 0 for quantum orders] else: @@ -219,9 +216,7 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, wrk_name = f"{workspace}_quantum_event_{n + 1}" partial_wrk_names.append(wrk_name) - self._fill_s_2d_workspace( - s_points=s_points[n], workspace=wrk_name, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_2d_workspace(s_points=s_points[n], workspace=wrk_name, species=species) GroupWorkspaces(InputWorkspaces=partial_wrk_names, OutputWorkspace=workspace) @@ -232,14 +227,12 @@ def _create_dummy_workspace(self, name): AnalysisDataService.addOrReplace(name, wrk) return wrk - def _fill_s_2d_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None): + def _fill_s_2d_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None): from mantid.api import NumericAxis from abins.constants import MILLI_EV_TO_WAVENUMBER - if protons_number is not None: - s_points = s_points * self.get_cross_section( - scattering=self._scale_by_cross_section, protons_number=protons_number, nucleons_number=nucleons_number - ) + if species is not None: + s_points = s_points * self.get_cross_section(scattering=self._scale_by_cross_section, species=species) n_q_values, n_freq_bins = s_points.shape n_q_bins = self._q_bins.size diff --git a/docs/source/release/v6.12.0/Indirect/Algorithms/Bugfixes/38448.rst b/docs/source/release/v6.12.0/Indirect/Algorithms/Bugfixes/38448.rst new file mode 100644 index 000000000000..3a07849a6847 --- /dev/null +++ b/docs/source/release/v6.12.0/Indirect/Algorithms/Bugfixes/38448.rst @@ -0,0 +1,13 @@ +The identification of isotopes has been improved in Abins/Abins2D. + +- Previously, a species would only be identified as the standard + isotopic mixture if the mass is very close to the Mantid reference + data. In some cases the values used by external phonon calculators + are significantly different and this could lead to misassignment. + (e.g. Zn in CASTEP.) Now, Abins will initially choose the _nearest_ + mass option between an isotope and the standard isotopic mixture. + +- Many isotopes lack cross-section data and would lead to NaN + intensities. Now, if a NaN cross-section is identified Abins + will either use the standard mixture data (if the mass is within + 0.01) or raise an error. diff --git a/scripts/abins/abins2.py b/scripts/abins/abins2.py index 4b0373ca99bd..6bea87526fa0 100644 --- a/scripts/abins/abins2.py +++ b/scripts/abins/abins2.py @@ -7,13 +7,15 @@ from typing import Dict, Optional +import numpy as np + from mantid.api import AlgorithmFactory, PythonAlgorithm, Progress from mantid.api import WorkspaceFactory, AnalysisDataService # noinspection PyProtectedMember from mantid.simpleapi import ConvertUnits, GroupWorkspaces import abins -from abins.abinsalgorithm import AbinsAlgorithm +from abins.abinsalgorithm import AbinsAlgorithm, AtomInfo from abins.logging import get_logger, Logger @@ -170,14 +172,13 @@ def PyExec(self): self.setProperty("OutputWorkspace", self._out_ws_name) prog_reporter.report("Group workspace with all required dynamical structure factors has been constructed.") - def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None): + def _fill_s_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None) -> None: """ Puts S into workspace(s). :param s_points: dynamical factor for the given atom - :param workspace: workspace to be filled with S - :param protons_number: number of protons in the given type fo atom - :param nucleons_number: number of nucleons in the given type of atom + :param workspace: workspace to be filled with S + :param species: Element/isotope data """ from abins.constants import FUNDAMENTALS, ONE_DIMENSIONAL_INSTRUMENTS, ONE_DIMENSIONAL_SPECTRUM @@ -187,15 +188,11 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, # only FUNDAMENTALS [data is 2d with one row] if s_points.shape[0] == FUNDAMENTALS: - self._fill_s_1d_workspace( - s_points=s_points[0], workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_1d_workspace(s_points=s_points[0], workspace=workspace, species=species) # total workspaces [data is 1d vector] elif len(s_points.shape) == ONE_DIMENSIONAL_SPECTRUM: - self._fill_s_1d_workspace( - s_points=s_points, workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_1d_workspace(s_points=s_points, workspace=workspace, species=species) # quantum order events (fundamentals or overtones + combinations for the given order) # [data is 2d table of S with a row for each quantum order] @@ -207,13 +204,11 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, wrk_name = f"{workspace}_quantum_event_{n + 1}" partial_wrk_names.append(wrk_name) - self._fill_s_1d_workspace( - s_points=s_points[n], workspace=wrk_name, protons_number=protons_number, nucleons_number=nucleons_number - ) + self._fill_s_1d_workspace(s_points=s_points[n], workspace=wrk_name, species=species) GroupWorkspaces(InputWorkspaces=partial_wrk_names, OutputWorkspace=workspace) - def _fill_s_1d_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None): + def _fill_s_1d_workspace(self, s_points=None, workspace=None, species: AtomInfo | None = None): """ Puts 1D S into workspace. :param protons_number: number of protons in the given type of atom @@ -221,10 +216,8 @@ def _fill_s_1d_workspace(self, s_points=None, workspace=None, protons_number=Non :param s_points: dynamical factor for the given atom :param workspace: workspace to be filled with S """ - if protons_number is not None: - s_points = s_points * self.get_cross_section( - scattering=self._scale_by_cross_section, protons_number=protons_number, nucleons_number=nucleons_number - ) + if species is not None: + s_points = s_points * self.get_cross_section(scattering=self._scale_by_cross_section, species=species) dim = 1 length = s_points.size diff --git a/scripts/abins/abinsalgorithm.py b/scripts/abins/abinsalgorithm.py index 1a501bbf0059..7cf30054b6cc 100644 --- a/scripts/abins/abinsalgorithm.py +++ b/scripts/abins/abinsalgorithm.py @@ -7,10 +7,13 @@ # Supporting functions for the Abins Algorithm that don't belong in # another part of AbinsModules. +import dataclasses +from functools import cached_property +from math import isnan import os from pathlib import Path import re -from typing import Dict, Iterable, List, Optional, Tuple, Union +from typing import Dict, Iterable, List, Literal, Tuple, Union import yaml @@ -24,10 +27,66 @@ from mantid.kernel import Atom, Direction, StringListValidator, StringArrayProperty, logger from mantid.simpleapi import CloneWorkspace, SaveAscii, Scale -import abins from abins.constants import AB_INITIO_FILE_EXTENSIONS, ALL_INSTRUMENTS, ATOM_PREFIX from abins.input.jsonloader import abins_supported_json_formats, JSONLoader from abins.instruments import get_instrument, Instrument +import abins.parameters + + +@dataclasses.dataclass +class AtomInfo: + symbol: str + mass: float + + @cached_property + def name(self): + if self.nucleons_number: + return f"{self.nucleons_number}{self.symbol}" + return self.symbol + + @property + def z_number(self): + return self._mantid_atom.z_number + + @property + def nucleons_number(self): + return self._mantid_atom.a_number + + @cached_property + def neutron_data(self): + return self._mantid_atom.neutron() + + @cached_property + def _mantid_atom(self): + from abins.constants import MASS_EPS + + nearest_int = int(round(self.mass)) + nearest_isotope = Atom(symbol=self.symbol, a_number=nearest_int) + standard_mix = Atom(symbol=self.symbol) + + if abs(nearest_isotope.mass - standard_mix.mass) < 1e-12: + # items are the same: standard mix is more likely to contain data + # (e.g. Atom('F', 19) has no neutron data but Atom('F') does) + return standard_mix + + if abs(self.mass - standard_mix.mass) < abs(self.mass - nearest_isotope.mass): + # Standard isotopic mixture, data should be available + return standard_mix + + if isnan(nearest_isotope.neutron().get("coh_scatt_length_real")): + logger.warning( + f"Nearest isotope to atomic mass {self.mass} is " + f"{nearest_int}{self.symbol}{nearest_isotope.z_number}, but " + f"neutron-scattering data is not available." + ) + + if abs(self.mass - standard_mix.mass) > MASS_EPS: + raise ValueError(f"Could not find suitable isotope data for {self.symbol} with mass {self.mass}.") + + logger.warning(f"Standard mass {standard_mix.mass} is close enough, will " f"use data for this isotope mixture.") + return standard_mix + + return nearest_isotope class AbinsAlgorithm: @@ -397,7 +456,7 @@ def create_workspaces(self, atoms_symbols=None, atom_numbers=None, *, s_data, at :returns: workspaces for list of atoms types, S for the particular type of atom """ - from abins.constants import FLOAT_TYPE, MASS_EPS, ONLY_ONE_MASS + from abins.constants import FLOAT_TYPE # Create appropriately-shaped arrays to be used in-place by _atom_type_s - avoid repeated slow instantiation shape = [max_quantum_order] @@ -412,7 +471,6 @@ def create_workspaces(self, atoms_symbols=None, atom_numbers=None, *, s_data, at if atoms_symbols is not None: for symbol in atoms_symbols: - sub = len(masses[symbol]) > ONLY_ONE_MASS or abs(Atom(symbol=symbol).mass - masses[symbol][0]) > MASS_EPS for m in masses[symbol]: result.extend( self._atom_type_s( @@ -423,7 +481,6 @@ def create_workspaces(self, atoms_symbols=None, atom_numbers=None, *, s_data, at element_symbol=symbol, temp_s_atom_data=temp_s_atom_data, s_atom_data=s_atom_data, - substitution=sub, ) ) if atom_numbers is not None: @@ -431,22 +488,26 @@ def create_workspaces(self, atoms_symbols=None, atom_numbers=None, *, s_data, at result.extend(self._atom_number_s(atom_number=atom_number, s_data=s_data, s_atom_data=s_atom_data, atoms_data=atoms_data)) return result - def _create_workspace(self, atom_name=None, s_points=None, optional_name="", protons_number=None, nucleons_number=None): + def _create_workspace(self, *, species: AtomInfo, s_points: np.ndarray, label: str | None = None): """ Creates workspace for the given frequencies and s_points with S data. After workspace is created it is rebined, scaled by cross-section factor and optionally multiplied by the user defined scaling factor. - :param atom_name: symbol of atom for which workspace should be created + :param species: Object identifying isotope (or mixture) :param s_points: S(Q, omega) - :param optional_name: optional part of workspace name + :param label: species-specific part of workspace name. If None, take from species object. :returns: workspace for the given frequency and S data - :param protons_number: number of protons in the given type fo atom - :param nucleons_number: number of nucleons in the given type of atom """ - ws_name = self._out_ws_name + "_" + atom_name + optional_name - self._fill_s_workspace(s_points=s_points, workspace=ws_name, protons_number=protons_number, nucleons_number=nucleons_number) + label = species.name if label is None else label + + ws_name = self._out_ws_name + "_" + label + self._fill_s_workspace( + species=species, + s_points=s_points, + workspace=ws_name, + ) return ws_name def _atom_number_s(self, *, atom_number, s_data, s_atom_data, atoms_data): @@ -474,8 +535,8 @@ def _atom_number_s(self, *, atom_number, s_data, s_atom_data, atoms_data): atom_workspaces = [] s_atom_data.fill(0.0) output_atom_label = "%s_%d" % (ATOM_PREFIX, atom_number) - symbol = atoms_data[atom_number - 1]["symbol"] - z_number = Atom(symbol=symbol).z_number + atom_data = atoms_data[atom_number - 1] + species = AtomInfo(symbol=atom_data["symbol"], mass=atom_data["mass"]) for i, order in enumerate(range(FUNDAMENTALS, self._max_event_order + 1)): s_atom_data[i] = s_data[atom_number - 1]["order_%s" % order] @@ -485,10 +546,12 @@ def _atom_number_s(self, *, atom_number, s_data, s_atom_data, atoms_data): atom_workspaces = [] atom_workspaces.append( self._create_workspace( - atom_name=output_atom_label, s_points=np.copy(total_s_atom_data), optional_name="_total", protons_number=z_number + species=species, + s_points=np.copy(total_s_atom_data), + label=output_atom_label + "_total", ) ) - atom_workspaces.append(self._create_workspace(atom_name=output_atom_label, s_points=np.copy(s_atom_data), protons_number=z_number)) + atom_workspaces.append(self._create_workspace(species=species, s_points=np.copy(s_atom_data), label=output_atom_label)) return atom_workspaces def _atom_type_s( @@ -500,7 +563,6 @@ def _atom_type_s( element_symbol=None, temp_s_atom_data=None, s_atom_data=None, - substitution=None, ): """ Helper function for calculating S for the given type of atom @@ -515,17 +577,17 @@ def _atom_type_s( information but is used in-place to save on time instantiating large arrays. :param s_atom_data: helper array to accumulate S (outer loop over atoms); does not transport information but is used in-place to save on time instantiating large arrays. - :param substitution: True if isotope substitution and False otherwise """ - from abins.constants import MASS_EPS + from abins.constants import FINE_MASS_EPS atom_workspaces = [] s_atom_data.fill(0.0) - element = Atom(symbol=element_symbol) + species = AtomInfo(symbol=element_symbol, mass=mass) for atom_index in range(num_atoms): - if atoms_data[atom_index]["symbol"] == element_symbol and abs(atoms_data[atom_index]["mass"] - mass) < MASS_EPS: + # TODO this mass_eps could be smaller, we are checking against known masses from data + if atoms_data[atom_index]["symbol"] == element_symbol and abs(atoms_data[atom_index]["mass"] - mass) < FINE_MASS_EPS: temp_s_atom_data.fill(0.0) for order in range(1, self._max_event_order + 1): @@ -537,35 +599,19 @@ def _atom_type_s( total_s_atom_data = np.sum(s_atom_data, axis=0) - nucleons_number = int(round(mass)) - - if substitution: - atom_workspaces.append( - self._create_workspace( - atom_name=str(nucleons_number) + element_symbol, - s_points=np.copy(total_s_atom_data), - optional_name="_total", - protons_number=element.z_number, - nucleons_number=nucleons_number, - ) - ) - atom_workspaces.append( - self._create_workspace( - atom_name=str(nucleons_number) + element_symbol, - s_points=np.copy(s_atom_data), - protons_number=element.z_number, - nucleons_number=nucleons_number, - ) - ) - else: - atom_workspaces.append( - self._create_workspace( - atom_name=element_symbol, s_points=np.copy(total_s_atom_data), optional_name="_total", protons_number=element.z_number - ) + atom_workspaces.append( + self._create_workspace( + species=species, + s_points=np.copy(total_s_atom_data), + label=f"{species.name}_total", ) - atom_workspaces.append( - self._create_workspace(atom_name=element_symbol, s_points=np.copy(s_atom_data), protons_number=element.z_number) + ) + atom_workspaces.append( + self._create_workspace( + species=species, + s_points=np.copy(s_atom_data), ) + ) return atom_workspaces @@ -616,7 +662,7 @@ def _create_total_workspace(self, partial_workspaces=None): s_atoms[:, i] += mtd[partial_ws].dataY(i) # create workspace with S - self._fill_s_workspace(s_atoms, total_workspace) + self._fill_s_workspace(s_points=s_atoms, workspace=total_workspace) # # Otherwise just repackage the workspace we have as the total else: @@ -642,26 +688,20 @@ def write_workspaces_to_ascii(scale: float = 1.0, *, ws_name: str) -> None: ) @staticmethod - def get_cross_section(scattering: str = "Total", nucleons_number: Optional[int] = None, *, protons_number: int) -> float: + def get_cross_section(scattering: Literal["Total", "Incoherent", "Coherent"], species: AtomInfo) -> float: """ Calculates cross section for the given element. :param scattering: Type of cross-section: 'Incoherent', 'Coherent' or 'Total' - :param protons_number: number of protons in the given type fo atom - :param nucleons_number: number of nucleons in the given type of atom + :param species: Data for atom/isotope type :returns: cross section for that element """ - if nucleons_number is not None: - try: - atom = Atom(a_number=nucleons_number, z_number=protons_number) - # isotopes are not implemented for all elements so use different constructor in that cases - except RuntimeError: - logger.warning(f"Could not find data for isotope {nucleons_number}, " f"using default values for {protons_number} protons.") - atom = Atom(z_number=protons_number) - else: - atom = Atom(z_number=protons_number) - scattering_keys = {"Incoherent": "inc_scatt_xs", "Coherent": "coh_scatt_xs", "Total": "tot_scatt_xs"} - return atom.neutron()[scattering_keys[scattering]] + cross_section = species.neutron_data[scattering_keys[scattering]] + + if isnan(cross_section): + raise ValueError(f"Found NaN cross-section for {species.symbol} with {species.nucleons_number} nucleons.") + + return cross_section @staticmethod def set_workspace_units(wrk, layout="1D", energy_units="cm-1"): diff --git a/scripts/abins/constants.py b/scripts/abins/constants.py index 4023fa6d3d3e..9049f3094dbe 100644 --- a/scripts/abins/constants.py +++ b/scripts/abins/constants.py @@ -317,10 +317,13 @@ ROTATIONS_AND_TRANSLATIONS = 6 -# This constant is used to check whether for the given atom mass averaged over all isotopes or mass of the -# specific isotope is used. +# This constant is used to check whether the standard atomic mass is "close enough", +# if closer isotope exists but lacks data MASS_EPS = 1e-2 # in amu units. +# Tolerance when distinguishing between mass types in the input data +FINE_MASS_EPS = 1e-6 + # this constant is used to check if in a system for the given symbol of an element all atoms with this symbol have # the same mass ONLY_ONE_MASS = 1 diff --git a/scripts/test/Abins/AbinsAlgorithmTest.py b/scripts/test/Abins/AbinsAlgorithmTest.py new file mode 100644 index 000000000000..729c9dd665bc --- /dev/null +++ b/scripts/test/Abins/AbinsAlgorithmTest.py @@ -0,0 +1,58 @@ +# Mantid Repository : https://github.com/mantidproject/mantid +# +# Copyright © 2024 ISIS Rutherford Appleton Laboratory UKRI, +# NScD Oak Ridge National Laboratory, European Spallation Source, +# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS +# SPDX - License - Identifier: GPL - 3.0 + + +import unittest + +# Import mantid.simpleapi first, otherwise we get circular import +import mantid.simpleapi # noqa: F401 + +from abins.abinsalgorithm import AbinsAlgorithm, AtomInfo + + +class AtomInfoTest(unittest.TestCase): + """Test the AtomInfo class""" + + def test_atom_info(self): + for args, expected in [ + # Non-standard isotope + (("Zn", 67.0), {"nucleons_number": 67, "name": "67Zn", "z_number": 30}), + # Round to standard mix + (("Zn", 65.4), {"nucleons_number": 0, "name": "Zn", "z_number": 30}), + # Choose standard mix when isotope with same mass is available + (("F", 19.0), {"nucleons_number": 0, "name": "F", "z_number": 9}), + ]: + atom_info = AtomInfo(*args) + for attr, value in expected.items(): + self.assertEqual(getattr(atom_info, attr), value) + + def test_bad_atom_info(self): + """Test that an error is raised if cross section data is missing""" + + # Zn65 is unstable and has no recorded cross section values + species = AtomInfo(symbol="Zn", mass=65.0) + + with self.assertRaisesRegex(ValueError, "Could not find suitable isotope data for Zn with mass 65.0."): + species.neutron_data + + +class AtomsDataTest(unittest.TestCase): + """Test static methods on AbinsAlgorithm""" + + def test_cross_section(self): + """Get cross section from nucleus information""" + + for scattering, nucleons_number, symbol, expected in [ + ("Incoherent", 67, "Zn", 0.28), + ("Coherent", 0, "Zn", 4.054), + ("Total", 0, "H", 82.02), + ]: + xc = AbinsAlgorithm.get_cross_section( + scattering=scattering, + species=AtomInfo(mass=float(nucleons_number), symbol=symbol), + ) + + self.assertAlmostEqual(xc, expected) diff --git a/scripts/test/Abins/CMakeLists.txt b/scripts/test/Abins/CMakeLists.txt index 28d4d325f905..ba978021026a 100644 --- a/scripts/test/Abins/CMakeLists.txt +++ b/scripts/test/Abins/CMakeLists.txt @@ -1,5 +1,6 @@ # Tests for Abins INS simulation set(TEST_PY_FILES + AbinsAlgorithmTest.py AbinsAtomsDataTest.py AbinsBroadeningTest.py AbinsPowderCalculatorTest.py