Refactor routine for estimating h ratio

When h ratio is not provided and hydrogen is present in the sample,
the routine for estimating the h ratio to the lowest mass is triggered.
Had to update this routine with the changes from converting vesuvio analysis into
a mantid algorithm.

src/mvesuvio/analysis_reduction.py

@@ -9,7 +9,7 @@
                             CloneWorkspace, DeleteWorkspace, VesuvioCalculateGammaBackground, \
                             VesuvioCalculateMS, Scale, RenameWorkspace, Minus, CreateSampleShape, \
                             VesuvioThickness, Integration, Divide, Multiply, DeleteWorkspaces, \
-                            CreateWorkspace
+                            CreateWorkspace, CreateSampleWorkspace
 
 from mvesuvio.util.analysis_helpers import loadConstants, numericalThirdDerivative
 
@@ -142,7 +142,7 @@ def PyExec(self):
         self.run()
 
     def _setup(self):
-        self._name = self.getProperty("InputWorkspace").value.name()
+        self._name = self.getPropertyValue("InputWorkspace")
         self._ip_file = self.getProperty("InstrumentParametersFile").value
         self._number_of_iterations = self.getProperty("NumberOfIterations").value
         self._mask_spectra = self.getProperty("InvalidDetectors").value 
@@ -191,17 +191,6 @@ def _setup(self):
         self._fit_profiles_workspaces = {}
 
 
-    def calculate_h_ratio(self):
-
-        if not np.isclose(min(self._masses), 1, atol=0.1):    # Hydrogen not present
-            return None
-
-        # Hydrogen present 
-        intensities = self._mean_intensity_ratios
-        sorted_intensities = intensities[np.argsort(self._masses)]
-
-        return sorted_intensities[0] / sorted_intensities[1] 
-
 
     def _update_workspace_data(self):
 

src/mvesuvio/analysis_routines.py

@@ -1,10 +1,12 @@
 # from .analysis_reduction import iterativeFitForDataReduction
 from mantid.api import AnalysisDataService
-from mantid.simpleapi import CreateEmptyTableWorkspace, mtd
+from mantid.simpleapi import CreateEmptyTableWorkspace, mtd, RenameWorkspace
 from mantid.api import AlgorithmFactory, AlgorithmManager
 import numpy as np
 
-from mvesuvio.util.analysis_helpers import fix_profile_parameters, loadRawAndEmptyWsFromUserPath, cropAndMaskWorkspace, NeutronComptonProfile
+from mvesuvio.util.analysis_helpers import fix_profile_parameters,  \
+                            loadRawAndEmptyWsFromUserPath, cropAndMaskWorkspace, \
+                            NeutronComptonProfile, calculate_h_ratio
 from mvesuvio.analysis_reduction import AnalysisRoutine
 from tests.testhelpers.calibration.algorithms import create_algorithm
 
@@ -39,7 +41,7 @@ def _create_analysis_object_from_current_interface(IC, running_tests=False):
     profiles_table = create_profiles_table(cropedWs.name()+"_initial_parameters", profiles)
 
     kwargs = {
-        "InputWorkspace": cropedWs,
+        "InputWorkspace": cropedWs.name(),
         "InputProfiles": profiles_table.name(),
         "InstrumentParametersFile": str(IC.InstrParsPath),
         "HRatioToLowestMass": IC.HToMassIdxRatio,
@@ -133,12 +135,14 @@ def run_joint_algs(back_alg, front_alg):
     assert incoming_means_table is not None, "Means table from backward routine not correctly accessed."
     assert h_ratio is not None, "H ratio from backward routine not correctly accesssed."
 
-    receiving_profiles_table = front_alg.getProperty("InputProfiles").value
+    receiving_profiles_table = mtd[front_alg.getPropertyValue("InputProfiles")]
 
-    fix_profile_parameters(incoming_means_table, receiving_profiles_table, h_ratio)
+    fixed_profiles_table = fix_profile_parameters(incoming_means_table, receiving_profiles_table, h_ratio)
 
+    # Update original profiles table
+    RenameWorkspace(fixed_profiles_table, receiving_profiles_table.name())
     # Even if the name is the same, need to trigger update
-    front_alg.setProperty("InputProfiles", receiving_profiles_table.name())
+    front_alg.setPropertyValue("InputProfiles", receiving_profiles_table.name())
 
     front_alg.execute()
     return
@@ -164,23 +168,26 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC):
 
     table_h_ratios = createTableWSHRatios()
 
-    backRoutine = _create_analysis_object_from_current_interface(bckwdIC)
-    frontRoutine = _create_analysis_object_from_current_interface(fwdIC)
+    back_alg = _create_analysis_object_from_current_interface(bckwdIC)
+    front_alg = _create_analysis_object_from_current_interface(fwdIC)
+
+    front_alg.execute()
+
+    means_table = mtd[front_alg.getPropertyValue("OutputMeansTable")]
+    current_ratio = calculate_h_ratio(means_table) 
 
-    frontRoutine.execute()
-    current_ratio = frontRoutine.calculate_h_ratio()
     table_h_ratios.addRow([current_ratio])
     previous_ratio = np.nan 
 
     while not np.isclose(current_ratio, previous_ratio, rtol=0.01):
 
-        backRoutine._h_ratio = current_ratio
-        backRoutine.execute()
-        frontRoutine.set_initial_profiles_from(backRoutine)
-        frontRoutine.execute()
+        back_alg.setProperty("HRatioToLowestMass", current_ratio)
+        run_joint_algs(back_alg, front_alg)
 
         previous_ratio = current_ratio
-        current_ratio = frontRoutine.calculate_h_ratio()
+
+        means_table = mtd[front_alg.getPropertyValue("OutputMeansTable")]
+        current_ratio = calculate_h_ratio(means_table) 
 
         table_h_ratios.addRow([current_ratio])
 
@@ -192,22 +199,9 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC):
 
 def createTableWSHRatios():
     table = CreateEmptyTableWorkspace(
-        OutputWorkspace="H_Ratios_Estimates"
+        OutputWorkspace="hydrogen_intensity_ratios_estimates"
     )
-    table.addColumn(type="float", name="H Ratio to lowest mass at each iteration")
+    table.addColumn(type="float", name="Hydrogen intensity ratio to lowest mass at each iteration")
     return table
 
 
-def isHPresent(masses) -> bool:
-    Hmask = np.abs(masses - 1) / 1 < 0.1  # H mass whithin 10% of 1 au
-
-    if np.any(Hmask):  # H present
-        print("\nH mass detected.\n")
-        assert (
-            len(Hmask) > 1
-        ), "When H is only mass present, run independent forward procedure, not joint."
-        assert Hmask[0], "H mass needs to be the first mass in masses and initPars."
-        assert sum(Hmask) == 1, "More than one mass very close to H were detected."
-        return True
-    else:
-        return False

src/mvesuvio/run_routine.py

@@ -8,11 +8,10 @@
     runIndependentIterativeProcedure,
     runJointBackAndForwardProcedure,
     runPreProcToEstHRatio,
-    createTableWSHRatios,
-    isHPresent,
 )
 
 from mantid.api import mtd
+import numpy as np
 
 def runRoutine(
     userCtr,
@@ -115,3 +114,18 @@ def checkInputs(crtIC):
 
     if (crtIC.procedure != "JOINT") & (crtIC.fitInYSpace is not None):
         assert crtIC.procedure == crtIC.fitInYSpace
+
+
+def isHPresent(masses) -> bool:
+    Hmask = np.abs(masses - 1) / 1 < 0.1  # H mass whithin 10% of 1 au
+
+    if np.any(Hmask):  # H present
+        print("\nH mass detected.\n")
+        assert (
+            len(Hmask) > 1
+        ), "When H is only mass present, run independent forward procedure, not joint."
+        assert Hmask[0], "H mass needs to be the first mass in masses and initPars."
+        assert sum(Hmask) == 1, "More than one mass very close to H were detected."
+        return True
+    else:
+        return False

src/mvesuvio/util/analysis_helpers.py

@@ -206,8 +206,7 @@ def fix_profile_parameters(incoming_means_table, receiving_profiles_table, h_rat
         p['width_bounds'] = str([p['width'] , p['width']])
 
     result_profiles_table = _convert_dict_to_table(profiles_dict)
-    RenameWorkspace(result_profiles_table, receiving_profiles_table.name())
-    return mtd[result_profiles_table.name()]
+    return result_profiles_table
 
 
 def _convert_table_to_dict(table):
@@ -234,3 +233,18 @@ def _get_lightest_profile(p_dict):
     profiles = [p for p in p_dict.values()]
     masses = [p['mass'] for p in p_dict.values()]
     return profiles[np.argmin(masses)]
+
+
+def calculate_h_ratio(means_table):
+
+    masses = means_table.column("mass")
+    intensities = np.array(means_table.column("mean_intensity"))
+
+    if not np.isclose(min(masses), 1, atol=0.1):    # Hydrogen not present
+        return None
+
+    # Hydrogen present 
+    sorted_intensities = intensities[np.argsort(masses)]
+
+    return sorted_intensities[0] / sorted_intensities[1] 
+

tests/unit/analysis/test_analysis_helpers.py

@@ -2,7 +2,8 @@
 import numpy as np
 import numpy.testing as nptest
 from mock import MagicMock
-from mvesuvio.util.analysis_helpers import extractWS, _convert_dict_to_table, fix_profile_parameters
+from mvesuvio.util.analysis_helpers import extractWS, _convert_dict_to_table,  \
+    fix_profile_parameters, calculate_h_ratio
 from mantid.simpleapi import CreateWorkspace, DeleteWorkspace
 
 
@@ -105,5 +106,13 @@ def test_fix_profile_parameters_without_H(self):
             {'label': '27.0', 'mass': 27.0, 'intensity': 0.3056112229824066, 'intensity_bounds': '[0, None]', 'width': 15.397000312805176, 'width_bounds': '[15.397, 15.397]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
         )
 
+
+    def test_calculate_h_ratio(self):
+        means_table_mock = MagicMock()
+        means_table_mock.column.side_effect = lambda x: [16, 1, 12] if x is "mass" else [0.1, 0.85, 0.05]
+        h_ratio = calculate_h_ratio(means_table_mock)
+        self.assertEqual(h_ratio, 0.85 / 0.05)
+
+
 if __name__ == "__main__":
     unittest.main()