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 404b22d8..ead6a60f 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
@@ -239,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
     )
 
 
@@ -453,6 +461,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 +617,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 +669,10 @@ def _neutron_compton_profiles(self, pars):
             / deltaQ
             * 0.72
         )
-        return intensities * (JOfY + FSE) * E0 * E0 ** (-0.92) * masses / deltaQ
+        scaling_factor = intensities * E0 * E0 ** (-0.92) * masses / deltaQ
+        JOfY *= scaling_factor
+        FSE *= scaling_factor
+        return JOfY+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
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