From b77fc39f88523d8ac70721e69bcf1e643c5949e7 Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Thu, 14 Nov 2024 15:54:58 +0000 Subject: [PATCH 1/4] Create workspaces for fse corrections Added workspaces to keep track of the FSE that were fitted to each spectra. These workspaces can then offer useful information to the scientist and can be optionally subtracted during the reduction procedure for the fit in y space. --- src/mvesuvio/analysis_reduction.py | 37 +++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 11 deletions(-) diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py index 404b22d8..85d11261 100644 --- a/src/mvesuvio/analysis_reduction.py +++ b/src/mvesuvio/analysis_reduction.py @@ -211,6 +211,12 @@ def _update_workspace_data(self): self._fit_profiles_workspaces[element] = self._create_emtpy_ncp_workspace(f'_{element}_ncp') self._fit_profiles_workspaces['total'] = self._create_emtpy_ncp_workspace(f'_total_ncp') + # Initialise workspaces for fitted ncp + self._fit_fse_workspaces = {} + for element in self._profiles_table.column("label"): + self._fit_fse_workspaces[element] = self._create_emtpy_ncp_workspace(f'_{element}_fse') + self._fit_fse_workspaces['total'] = self._create_emtpy_ncp_workspace(f'_total_fse') + # Initialise empty means self._mean_widths = None self._std_widths = None @@ -453,6 +459,12 @@ def _create_summed_workspaces(self): OutputWorkspace=ws.name() + "_sum" ) + for ws in self._fit_fse_workspaces.values(): + SumSpectra( + InputWorkspace=ws.name(), + OutputWorkspace=ws.name() + "_sum" + ) + def _set_means_and_std(self): """Calculate mean widths and intensities from tableWorkspace""" @@ -603,30 +615,32 @@ def _fit_neutron_compton_profiles_to_row(self): self.log().notice(' '.join(str(tableRow).split(",")).replace('[', '').replace(']', '')) # Pass fit profiles into workspaces - individual_ncps = self._neutron_compton_profiles(fitPars) - for ncp, element in zip(individual_ncps, self._profiles_table.column("label")): + ncp_for_each_mass, fse_for_each_mass = self._neutron_compton_profiles(fitPars) + for ncp, fse, element in zip(ncp_for_each_mass, fse_for_each_mass, self._profiles_table.column("label")): self._fit_profiles_workspaces[element].dataY(self._row_being_fit)[:] = ncp + self._fit_fse_workspaces[element].dataY(self._row_being_fit)[:] = fse - self._fit_profiles_workspaces['total'].dataY(self._row_being_fit)[:] = np.sum(individual_ncps, axis=0) + self._fit_profiles_workspaces['total'].dataY(self._row_being_fit)[:] = np.sum(ncp_for_each_mass, axis=0) + self._fit_fse_workspaces['total'].dataY(self._row_being_fit)[:] = np.sum(fse_for_each_mass, axis=0) return def errorFunction(self, pars): """Error function to be minimized, in TOF space""" - ncpForEachMass = self._neutron_compton_profiles(pars) - ncpTotal = np.sum(ncpForEachMass, axis=0) + ncp_for_each_mass, fse_for_each_mass = self._neutron_compton_profiles(pars) + ncp_total = np.sum(ncp_for_each_mass, axis=0) # Ignore any masked values from Jackknife or masked tof range zerosMask = self._dataY[self._row_being_fit] == 0 - ncpTotal = ncpTotal[~zerosMask] - dataY = self._dataY[self._row_being_fit, ~zerosMask] - dataE = self._dataE[self._row_being_fit, ~zerosMask] + ncp_total = ncp_total[~zerosMask] + data_y = self._dataY[self._row_being_fit, ~zerosMask] + data_e = self._dataE[self._row_being_fit, ~zerosMask] if np.all(self._dataE[self._row_being_fit] == 0): # When errors not present - return np.sum((ncpTotal - dataY) ** 2) + return np.sum((ncp_total - data_y) ** 2) - return np.sum((ncpTotal - dataY) ** 2 / dataE**2) + return np.sum((ncp_total - data_y) ** 2 / data_e**2) def _neutron_compton_profiles(self, pars): @@ -653,7 +667,8 @@ def _neutron_compton_profiles(self, pars): / deltaQ * 0.72 ) - return intensities * (JOfY + FSE) * E0 * E0 ** (-0.92) * masses / deltaQ + NCP = intensities * (JOfY + FSE) * E0 * E0 ** (-0.92) * masses / deltaQ + return NCP, FSE def caculateResolutionForEachMass(self, centers): From 269cc17cc6431b7a5441a5b273d2186483b6603b Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Fri, 15 Nov 2024 15:47:41 +0000 Subject: [PATCH 2/4] Add option for subtracting FSE --- src/mvesuvio/analysis_fitting.py | 21 +++++++++++++++++++++ src/mvesuvio/analysis_reduction.py | 10 +++++++--- src/mvesuvio/config/analysis_inputs.py | 3 ++- 3 files changed, 30 insertions(+), 4 deletions(-) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index fd1867c0..3d5e4cb1 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -8,6 +8,7 @@ from iminuit.util import make_func_code, describe import jacobi import time +import re from mvesuvio.util import handle_config @@ -27,6 +28,10 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): wsTOFMass0 = subtractAllMassesExceptFirst(IC, wsTOF, ncpForEachMass) + if yFitIC.subtractFSE: + wsFSEMass0 = find_ws_name_fse_first_mass(IC) + wsTOFMass0 = Minus(wsTOFMass0, wsFSEMass0, OutputWorkspace=wsTOFMass0.name() + "_fse") + wsJoY, wsJoYAvg = ySpaceReduction( wsTOFMass0, IC.masses[0], yFitIC, ncpForEachMass[:, 0, :] ) @@ -48,6 +53,22 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): return yfitResults +def find_ws_name_fse_first_mass(ic): + # Some dirty implementation to be deleted in the future + ws_names_fse = [] + ws_masses = [] + prefix = ic.name+'_'+str(ic.noOfMSIterations) + for ws_name in mtd.getObjectNames(): + if ws_name.startswith(prefix) and ws_name.endswith('fse'): + name_ending = ws_name.replace(prefix, "") + match = re.search(r'\d+\.?\d*', name_ending) + if match: # If float found in ws name ending + ws_masses.append(float(match.group())) + ws_names_fse.append(ws_name) + + return ws_names_fse[np.argmin(ws_masses)] + + def extractNCPFromWorkspaces(wsFinal, ic): """Extra function to extract ncps from loaded ws in mantid.""" diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py index 85d11261..5eeace81 100644 --- a/src/mvesuvio/analysis_reduction.py +++ b/src/mvesuvio/analysis_reduction.py @@ -245,8 +245,10 @@ def _create_emtpy_ncp_workspace(self, suffix): DataY=np.zeros(self._dataY.size), DataE=np.zeros(self._dataE.size), Nspec=self._workspace_being_fit.getNumberHistograms(), + UnitX="TOF", # I had hoped for a method like .XUnit() but alas OutputWorkspace=self._workspace_being_fit.name()+suffix, - ParentWorkspace=self._workspace_being_fit + ParentWorkspace=self._workspace_being_fit, + Distribution=True ) @@ -667,8 +669,10 @@ def _neutron_compton_profiles(self, pars): / deltaQ * 0.72 ) - NCP = intensities * (JOfY + FSE) * E0 * E0 ** (-0.92) * masses / deltaQ - return NCP, FSE + scaling_factor = intensities * E0 * E0 ** (-0.92) * masses / deltaQ + JoY *= scaling_factor + FSE *= scaling_factor + return JoY+FSE, FSE def caculateResolutionForEachMass(self, centers): diff --git a/src/mvesuvio/config/analysis_inputs.py b/src/mvesuvio/config/analysis_inputs.py index 448986a7..af287c78 100644 --- a/src/mvesuvio/config/analysis_inputs.py +++ b/src/mvesuvio/config/analysis_inputs.py @@ -86,7 +86,8 @@ class ForwardAnalysisInputs(SampleParameters): @dataclass class YSpaceFitInputs: showPlots = True - symmetrisationFlag = True + symmetrisationFlag = False + subtractFSE = True rebinParametersForYSpaceFit = "-25, 0.5, 25" # Needs to be symetric fitModel = "SINGLE_GAUSSIAN" # Options: 'SINGLE_GAUSSIAN', 'GC_C4', 'GC_C6', 'GC_C4_C6', 'DOUBLE_WELL', 'ANSIO_GAUSSIAN', 'Gaussian3D' runMinos = True From 0a523c84a434a0356b7fe5ef578d2e23643b490f Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Fri, 15 Nov 2024 15:52:43 +0000 Subject: [PATCH 3/4] Add subtract FSE option to test inputs --- tests/data/analysis/inputs/analysis_test.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/data/analysis/inputs/analysis_test.py b/tests/data/analysis/inputs/analysis_test.py index 99d31d13..dbf660ec 100644 --- a/tests/data/analysis/inputs/analysis_test.py +++ b/tests/data/analysis/inputs/analysis_test.py @@ -87,6 +87,7 @@ class ForwardAnalysisInputs(SampleParameters): class YSpaceFitInputs: showPlots = False symmetrisationFlag = True + subtractFSE = False rebinParametersForYSpaceFit = "-20, 0.5, 20" # Needs to be symetric fitModel = "SINGLE_GAUSSIAN" runMinos = True From 2b2334848add964cc1aaaf47913b8cb71465e7de Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Fri, 15 Nov 2024 16:09:40 +0000 Subject: [PATCH 4/4] Fix syntax error --- src/mvesuvio/analysis_reduction.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py index 5eeace81..ead6a60f 100644 --- a/src/mvesuvio/analysis_reduction.py +++ b/src/mvesuvio/analysis_reduction.py @@ -670,9 +670,9 @@ def _neutron_compton_profiles(self, pars): * 0.72 ) scaling_factor = intensities * E0 * E0 ** (-0.92) * masses / deltaQ - JoY *= scaling_factor + JOfY *= scaling_factor FSE *= scaling_factor - return JoY+FSE, FSE + return JOfY+FSE, FSE def caculateResolutionForEachMass(self, centers):