From 787f3f0e9177506846c9837b201628a707f9d5f4 Mon Sep 17 00:00:00 2001
From: Jess Farmer <jessica.farmer@stfc.ac.uk>
Date: Thu, 27 Apr 2023 15:53:03 +0000
Subject: [PATCH] allow fitting of multiple peaks and add new system test
---
.../calibrate_vesuvio_analysis.py | 65 +++++++++++++------
.../calibrate_vesuvio_fit.py | 40 +++++++-----
.../tests/system/test_system_analysis.py | 26 +++++++-
.../tests/unit/test_calibrate_vesuvio_fit.py | 27 ++++----
4 files changed, 108 insertions(+), 50 deletions(-)
diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py
index e5431294..e9a28b79 100644
--- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py
+++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_analysis.py
@@ -136,7 +136,7 @@ def PyExec(self):
self._calculate_final_energy(E1_peak_fits_back, EVSGlobals.BACKSCATTERING_RANGE, self._shared_parameter_fit_type)
# calibrate L1 for backscattering detectors based on the averaged E1 value and calibrated theta values
- self._calculate_final_flight_path(E1_peak_fits_back[0], EVSGlobals.BACKSCATTERING_RANGE)
+ self._calculate_final_flight_path(E1_peak_fits_back, EVSGlobals.BACKSCATTERING_RANGE, self._shared_parameter_fit_type)
# calibrate L1 for frontscattering detectors based on the averaged E1 value and calibrated theta values
E1_fit_front = self._current_workspace + '_E1_front'
@@ -146,7 +146,7 @@ def PyExec(self):
PeakType='Recoil', SharedParameterFitType=self._shared_parameter_fit_type)
E1_peak_fits_front = mtd[self._current_workspace + '_E1_front_Peak_Parameters'].getNames()
- self._calculate_final_flight_path(E1_peak_fits_front[0], EVSGlobals.FRONTSCATTERING_RANGE)
+ self._calculate_final_flight_path(E1_peak_fits_front, EVSGlobals.FRONTSCATTERING_RANGE, self._shared_parameter_fit_type)
# make the fitted parameters for this iteration the input to the next iteration.
table_group.append(self._current_workspace)
@@ -234,39 +234,66 @@ def _calculate_incident_flight_path(self, table_name, spec_list):
DeleteWorkspace(L0_param_table)
- def _calculate_final_flight_path(self, peak_table, spec_list):
+ def _calculate_final_flight_path(self, peak_table, spec_list, fit_type):
"""
Calculate the final flight path using the values for energy.
This also uses the old value for L1 loaded from the parameter file.
@param spec_list - spectrum range to calculate t0 for.
"""
+ if fit_type == 'Both':
+ indv_peak_table = peak_table[0]
+ glob_peak_table = peak_table[1]
+ E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list)
+ E1 *= EVSGlobals.MEV_CONVERSION
+ global_E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'global_E1', spec_list)
+ global_E1 *= EVSGlobals.MEV_CONVERSION
+ elif fit_type == 'Shared':
+ glob_peak_table = peak_table[0]
+ global_E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'global_E1', spec_list)
+ global_E1 *= EVSGlobals.MEV_CONVERSION
+ else:
+ indv_peak_table = peak_table[0]
+ E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list)
+ E1 *= EVSGlobals.MEV_CONVERSION
- E1 = EVSMiscFunctions.read_table_column(self._current_workspace, 'E1', spec_list)
t0 = EVSMiscFunctions.read_table_column(self._current_workspace, 't0', spec_list)
t0_error = EVSMiscFunctions.read_table_column(self._current_workspace, 't0_Err', spec_list)
L0 = EVSMiscFunctions.read_table_column(self._current_workspace, 'L0', spec_list)
theta = EVSMiscFunctions.read_table_column(self._current_workspace, 'theta', spec_list)
r_theta = EVSMiscFunctions.calculate_r_theta(self._sample_mass, theta)
- peak_centres = EVSMiscFunctions.read_fitting_result_table_column(peak_table, 'f1.LorentzPos', spec_list)
- peak_centres_errors = EVSMiscFunctions.read_fitting_result_table_column(peak_table, 'f1.LorentzPos_Err', spec_list)
- invalid_spectra = EVSMiscFunctions.identify_invalid_spectra(peak_table, peak_centres, peak_centres_errors, spec_list)
- peak_centres[invalid_spectra] = np.nan
+ if fit_type != 'Shared':
+ peak_centres = EVSMiscFunctions.read_fitting_result_table_column(indv_peak_table, 'f1.LorentzPos', spec_list)
+ peak_centres_errors = EVSMiscFunctions.read_fitting_result_table_column(indv_peak_table, 'f1.LorentzPos_Err', spec_list)
+ invalid_spectra = EVSMiscFunctions.identify_invalid_spectra(indv_peak_table, peak_centres, peak_centres_errors, spec_list)
+ peak_centres[invalid_spectra] = np.nan
+
+ print(f'Invalid Spectra Index Found and Marked NAN: {invalid_spectra.flatten()} from Spectra Index List:'
+ f'{[ x -3 for x in spec_list]}')
+
+ delta_t = (peak_centres - t0) / 1e+6
+ delta_t_error = t0_error / 1e+6
+ v1 = np.sqrt( 2 * E1 /scipy.constants.m_n)
+ L1 = v1 * delta_t - L0 * r_theta
+ L1_error = v1 * delta_t_error
- print(f'Invalid Spectra Index Found and Marked NAN: {invalid_spectra.flatten()} from Spectra Index List:'
- f'{[ x -3 for x in spec_list]}')
+ self._set_table_column(self._current_workspace, 'L1', L1, spec_list)
+ self._set_table_column(self._current_workspace, 'L1_Err', L1_error, spec_list)
- delta_t = (peak_centres - t0) / 1e+6
- delta_t_error = t0_error / 1e+6
+ if fit_type != 'Individual':
+ peak_centre = EVSMiscFunctions.read_fitting_result_table_column(glob_peak_table, 'f1.LorentzPos', [0])
+ peak_centre_list = np.empty(spec_list[1] + 1 -spec_list[0])
+ peak_centre_list.fill(float(peak_centre))
- E1 *= EVSGlobals.MEV_CONVERSION
- v1 = np.sqrt( 2 * E1 /scipy.constants.m_n)
- L1 = v1 * delta_t - L0 * r_theta
- L1_error = v1 * delta_t_error
+ delta_t = (peak_centre_list - t0) / 1e+6
+ delta_t_error = t0_error / 1e+6
+ global_v1 = np.sqrt( 2 * global_E1 /scipy.constants.m_n)
+ global_L1 = global_v1 * delta_t - L0 * r_theta
+ global_L1_error = global_v1 * delta_t_error
- self._set_table_column(self._current_workspace, 'L1', L1, spec_list)
- self._set_table_column(self._current_workspace, 'L1_Err', L1_error, spec_list)
+ self._set_table_column(self._current_workspace, 'global_L1', global_L1, spec_list)
+ self._set_table_column(self._current_workspace, 'global_L1_Err', global_L1_error, spec_list)
def _calculate_scattering_angle(self, table_name, spec_list):
"""
@@ -383,7 +410,7 @@ def _calculate_final_energy(self, peak_table, spec_list, fit_type):
self._set_table_column(self._current_workspace, 'global_E1', global_E1)
self._set_table_column(self._current_workspace, 'global_E1_Err', global_E1_error)
- else:
+ else:
mean_E1 = np.empty(spec_range)
E1_error = np.empty(spec_range)
diff --git a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py
index 8275cada..510e01a4 100644
--- a/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py
+++ b/unpackaged/vesuvio_calibration/calibration_scripts/calibrate_vesuvio_fit.py
@@ -546,13 +546,12 @@ def _fit_peaks(self):
"""
estimated_peak_positions_all_peaks = self._estimate_peak_positions()
+ num_estimated_peaks, num_spectra = estimated_peak_positions_all_peaks.shape
+ self._output_parameter_tables = []
+ self._peak_fit_workspaces = []
if self._shared_parameter_fit_type != "Shared":
- num_estimated_peaks, num_spectra = estimated_peak_positions_all_peaks.shape
-
self._prog_reporter = Progress(self, 0.0, 1.0, num_estimated_peaks*num_spectra)
- self._output_parameter_tables = []
- self._peak_fit_workspaces = []
for peak_index, estimated_peak_positions in enumerate(estimated_peak_positions_all_peaks):
self._peak_fit_workspaces_by_spec = []
@@ -570,18 +569,21 @@ def _fit_peaks(self):
GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters')
if self._shared_parameter_fit_type != "Individual":
- estimated_peak_position = np.mean(estimated_peak_positions_all_peaks)
- output_parameter_table_name = self._output_workspace_name + '_Shared_Peak_Parameters'
- output_parameter_table_headers = self._create_parameter_table_and_output_headers(output_parameter_table_name)
+ for peak_index, estimated_peak_position in enumerate(np.mean(estimated_peak_positions_all_peaks,1)):
- self._fit_shared_peak(self._spec_list[0], estimated_peak_position, output_parameter_table_name,
- output_parameter_table_headers)
+ output_parameter_table_name = self._output_workspace_name + '_Shared_Peak_%d_Parameters' % peak_index
+ output_parameter_table_headers = self._create_parameter_table_and_output_headers(output_parameter_table_name)
+
+ fit_workspace_name = self._fit_shared_peak(peak_index, self._spec_list[0], estimated_peak_position, output_parameter_table_name,
+ output_parameter_table_headers)
- if self._shared_parameter_fit_type == 'Both':
self._output_parameter_tables.append(output_parameter_table_name)
+ self._peak_fit_workspaces.append(fit_workspace_name)
+
+ if self._shared_parameter_fit_type == 'Both':
GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters')
else:
- GroupWorkspaces(output_parameter_table_name, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters')
+ GroupWorkspaces(self._output_parameter_tables, OutputWorkspace=self._output_workspace_name + '_Peak_Parameters')
def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_table_name,
output_parameter_table_headers):
@@ -605,8 +607,8 @@ def _fit_peak(self, peak_index, spec_index, peak_position, output_parameter_tabl
return fit_workspace_name
- def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_name, output_parameter_table_headers):
- peak_params = self._find_peaks_and_output_params(spec_range, peak_position)
+ def _fit_shared_peak(self, peak_index, spec_range, peak_position, output_parameter_table_name, output_parameter_table_headers):
+ peak_params = self._find_peaks_and_output_params(spec_range, peak_position, peak_index)
fit_func = self._build_function_string(peak_params)
start_str = 'composite=MultiDomainFunction, NumDeriv=1;'
validSpecs = mtd[self._sample]
@@ -621,7 +623,7 @@ def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_nam
ties = ','.join(f'f{i}.f1.{p}=f0.f1.{p}' for p in self._func_param_names.values() for i in range(1,n_valid_specs))
func_string = composite_func + f';ties=({ties})'
xmin, xmax = self._find_fit_x_window(peak_params)
- fit_output_name = '__' + self._output_workspace_name + '_Peak_0'
+ fit_output_name = '__' + self._output_workspace_name + '_Peak_%d' % peak_index
# create new workspace for each spectra
x = validSpecs.readX(0)
@@ -646,13 +648,17 @@ def _fit_shared_peak(self, spec_range, peak_position, output_parameter_table_nam
self._output_fit_params_to_table_ws(0, fit_params, output_parameter_table_name,
output_headers)
+
+ fit_workspace_name = fws.name()
self._del_fit_workspaces(ncm, fit_params, fws)
- def _find_peaks_and_output_params(self, spec_number, peak_position, peak_index=None, spec_index=None):
- if spec_index is not None and peak_index is not None:
+ return fit_workspace_name
+
+ def _find_peaks_and_output_params(self, spec_number, peak_position, peak_index, spec_index=None):
+ if spec_index is not None:
peak_table_name = '__' + self._sample + '_peaks_table_%d_%d' % (peak_index, spec_index)
else:
- peak_table_name = '__' + self._sample + '_peak_table'
+ peak_table_name = '__' + self._sample + '_peak_table_%d' % (peak_index)
find_peak_params = self._get_find_peak_parameters(spec_number, [peak_position])
FindPeaks(InputWorkspace=self._sample, WorkspaceIndex=spec_index, PeaksList=peak_table_name, **find_peak_params)
if mtd[peak_table_name].rowCount() == 0:
diff --git a/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py b/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py
index 59c0b88c..cf986d55 100644
--- a/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py
+++ b/unpackaged/vesuvio_calibration/tests/system/test_system_analysis.py
@@ -190,9 +190,27 @@ def test_copper_with_individual_and_global_fit(self, load_file_mock):
detector_specific_r_tols = {"L1": {116: 0.65, 170: 0.75}}
detector_specific_r_tols["L1"].update({k: TestConstants.INVALID_DETECTOR for k in [138, 141, 146, 147, 156, 158, 160, 163, 164,
165, 167, 168, 169, 170, 182, 191, 192]})
- self._assert_parameters_match_expected(params_table, detector_specific_r_tols)
+ self._assert_parameters_match_expected(params_table, detector_specific_r_tols, "Both")
self._assert_E1_parameters_match_expected(params_table, 3e-3, "Both")
+ @patch('unpackaged.vesuvio_calibration.calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._load_file')
+ def test_copper_with_global_fit(self, load_file_mock):
+ self._setup_copper_test()
+ self._output_workspace = "copper_analysis_test"
+
+ load_file_mock.side_effect = self._load_file_side_effect
+
+ self._create_evs_calibration_alg()
+ self._alg.setProperty("SharedParameterFitType", "Shared")
+ params_table = self._execute_evs_calibration_analysis()
+
+ # Specify detectors tolerances set by user, then update with those to mask as invalid.
+ detector_specific_r_tols = {"L1": {116: 0.65, 170: 0.75}}
+ detector_specific_r_tols["L1"].update({k: TestConstants.INVALID_DETECTOR for k in [138, 141, 146, 147, 156, 158, 160, 163, 164,
+ 165, 167, 168, 169, 170, 182, 191, 192]})
+ self._assert_parameters_match_expected(params_table, detector_specific_r_tols, "Shared")
+ self._assert_E1_parameters_match_expected(params_table, 3e-3, "Shared")
+
def _assert_theta_parameters_match_expected(self, params_table, rel_tolerance):
thetas = params_table.column('theta')
actual_thetas = self._calibrated_params['theta']
@@ -225,10 +243,12 @@ def _assert_E1_parameters_match_expected(self, params_table, rel_tolerance, fit_
global_E1 = params_table.column('global_E1')[0]
self.assertAlmostEqual(global_E1, EVSGlobals.ENERGY_ESTIMATE, delta=EVSGlobals.ENERGY_ESTIMATE*rel_tolerance)
- def _assert_parameters_match_expected(self, params_table, tolerances=None):
+ def _assert_parameters_match_expected(self, params_table, tolerances=None, fit_type="Individual"):
rel_tol_theta, rel_tol_L1 = self._extract_tolerances(tolerances)
theta_errors = self._assert_theta_parameters_match_expected(params_table, rel_tol_theta)
- L1_errors = self._assert_L1_parameters_match_expected(params_table, rel_tol_L1)
+ L1_errors = None
+ if fit_type != 'Shared':
+ L1_errors = self._assert_L1_parameters_match_expected(params_table, rel_tol_L1)
if theta_errors or L1_errors:
raise AssertionError(f"Theta: {theta_errors})\n L1: {L1_errors}")
diff --git a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py
index aaaa0b2c..ae79d63f 100644
--- a/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py
+++ b/unpackaged/vesuvio_calibration/tests/unit/test_calibrate_vesuvio_fit.py
@@ -660,29 +660,32 @@ def test_fit_peaks_individual(self, group_workspaces_mock):
group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters'],
OutputWorkspace='output_ws_name_Peak_Parameters')
- @patch('unpackaged.vesuvio_calibration.calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._fit_shared_peak')
@patch('unpackaged.vesuvio_calibration.calibration_scripts.calibrate_vesuvio_fit.GroupWorkspaces')
- def test_fit_peaks_shared(self, group_workspaces_mock, fit_shared_peak_mock):
+ def test_fit_peaks_shared(self, group_workspaces_mock):
alg = EVSCalibrationFit()
alg._estimate_peak_positions = MagicMock(return_value=np.asarray([[5, 10, 15], [2.5, 7.5, 10.5]]))
alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c'])
- alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}')
+ alg._fit_shared_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}')
alg._output_workspace_name = 'output_ws_name'
alg._shared_parameter_fit_type = 'Shared'
alg._spec_list = [3,4,5]
alg._fit_peaks()
- alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Shared_Peak_Parameters')])
- group_workspaces_mock.assert_called_once_with('output_ws_name_Shared_Peak_Parameters',
+ self.assertEqual(['fit_ws_0', 'fit_ws_1'],
+ alg._peak_fit_workspaces)
+ alg._create_parameter_table_and_output_headers.assert_has_calls([call('output_ws_name_Shared_Peak_0_Parameters'),
+ call('output_ws_name_Shared_Peak_1_Parameters')])
+ group_workspaces_mock.assert_called_once_with(['output_ws_name_Shared_Peak_0_Parameters','output_ws_name_Shared_Peak_1_Parameters'],
OutputWorkspace='output_ws_name_Peak_Parameters')
- fit_shared_peak_mock.assert_called_once_with(3, np.mean(alg._estimate_peak_positions()), 'output_ws_name_Shared_Peak_Parameters', ['a', 'b', 'c'])
+ alg._fit_shared_peak.assert_has_calls([call(0, 3, 10, 'output_ws_name_Shared_Peak_0_Parameters', ['a', 'b', 'c']),
+ call(1, 3, 41/6, 'output_ws_name_Shared_Peak_1_Parameters', ['a', 'b', 'c'])])
- @patch('unpackaged.vesuvio_calibration.calibration_scripts.calibrate_vesuvio_fit.EVSCalibrationFit._fit_shared_peak')
@patch('unpackaged.vesuvio_calibration.calibration_scripts.calibrate_vesuvio_fit.GroupWorkspaces')
- def test_fit_peaks_both(self, group_workspaces_mock, fit_shared_peak_mock):
+ def test_fit_peaks_both(self, group_workspaces_mock):
alg = EVSCalibrationFit()
alg._estimate_peak_positions = MagicMock(return_value=np.asarray([[5, 10, 15], [2.5, 7.5, 10.5]]))
alg._create_parameter_table_and_output_headers = MagicMock(return_value=['a', 'b', 'c'])
alg._fit_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}_{b}')
+ alg._fit_shared_peak = MagicMock(side_effect=lambda a, b, c, d, e: f'fit_ws_{a}')
alg._output_workspace_name = 'output_ws_name'
alg._shared_parameter_fit_type = 'Both'
alg._spec_list = [3,4,5]
@@ -695,11 +698,13 @@ def test_fit_peaks_both(self, group_workspaces_mock, fit_shared_peak_mock):
call(1, 0, 2.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']),
call(1, 1, 7.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c']),
call(1, 2, 10.5, 'output_ws_name_Peak_1_Parameters', ['a', 'b', 'c'])])
- self.assertEqual([['fit_ws_0_0', 'fit_ws_0_1', 'fit_ws_0_2'], ['fit_ws_1_0', 'fit_ws_1_1', 'fit_ws_1_2']],
+ self.assertEqual([['fit_ws_0_0', 'fit_ws_0_1', 'fit_ws_0_2'], ['fit_ws_1_0', 'fit_ws_1_1', 'fit_ws_1_2'], 'fit_ws_0', 'fit_ws_1'],
alg._peak_fit_workspaces)
group_workspaces_mock.assert_called_once_with(['output_ws_name_Peak_0_Parameters', 'output_ws_name_Peak_1_Parameters',
- 'output_ws_name_Shared_Peak_Parameters'], OutputWorkspace='output_ws_name_Peak_Parameters')
- fit_shared_peak_mock.assert_called_once_with(3, np.mean(alg._estimate_peak_positions()), 'output_ws_name_Shared_Peak_Parameters', ['a', 'b', 'c'])
+ 'output_ws_name_Shared_Peak_0_Parameters','output_ws_name_Shared_Peak_1_Parameters'],
+ OutputWorkspace='output_ws_name_Peak_Parameters')
+ alg._fit_shared_peak.assert_has_calls([call(0, 3, 10, 'output_ws_name_Shared_Peak_0_Parameters', ['a', 'b', 'c']),
+ call(1, 3, 41/6, 'output_ws_name_Shared_Peak_1_Parameters', ['a', 'b', 'c'])])
@staticmethod
def _setup_alg_mocks_fit_peak():