diff --git a/mvesuvio/config/analysis_inputs.py b/mvesuvio/config/analysis_inputs.py
index b005d7a6..dac2fd3e 100644
--- a/mvesuvio/config/analysis_inputs.py
+++ b/mvesuvio/config/analysis_inputs.py
@@ -139,7 +139,8 @@ def __init__(self, ipFilesPath):
 
 class YSpaceFitInitialConditions:
     showPlots = True
-    symmetrisationFlag = True
+    subtractFSE = True
+    symmetrisationFlag = False 
     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/mvesuvio/system_tests/test_config/expr_for_tests.py b/mvesuvio/system_tests/test_config/expr_for_tests.py
index 4f9dd111..9fb59ac0 100644
--- a/mvesuvio/system_tests/test_config/expr_for_tests.py
+++ b/mvesuvio/system_tests/test_config/expr_for_tests.py
@@ -108,6 +108,7 @@ def __init__(self, ipFilesPath):
 
 class YSpaceFitInitialConditions:
     showPlots = False
+    subtractFSE = False
     symmetrisationFlag = True
     rebinParametersForYSpaceFit = "-20, 0.5, 20"  # Needs to be symetric
     fitModel = "SINGLE_GAUSSIAN"
diff --git a/mvesuvio/vesuvio_analysis/analysis_functions.py b/mvesuvio/vesuvio_analysis/analysis_functions.py
index b56c5000..4cc882f0 100644
--- a/mvesuvio/vesuvio_analysis/analysis_functions.py
+++ b/mvesuvio/vesuvio_analysis/analysis_functions.py
@@ -249,7 +249,7 @@ def fitNcpToWorkspace(IC, ws):
     )
     createTableWSForFitPars(ws.name(), IC.noOfMasses, arrFitPars)
     arrBestFitPars = arrFitPars[:, 1:-2]
-    ncpForEachMass, ncpTotal = calculateNcpArr(
+    ncpForEachMass, ncpTotal, fseForEachMass = calculateNcpArr(
         IC,
         arrBestFitPars,
         resolutionPars,
@@ -258,6 +258,7 @@ def fitNcpToWorkspace(IC, ws):
         ySpacesForEachMass,
     )
     ncpSumWSs = createNcpWorkspaces(ncpForEachMass, ncpTotal, ws, IC)
+    createFseWorkspaces(fseForEachMass, ws, IC)
 
     wsDataSum = SumSpectra(InputWorkspace=ws, OutputWorkspace=ws.name() + "_Sum")
     plotSumNCPFits(wsDataSum, *ncpSumWSs, IC)
@@ -437,8 +438,9 @@ def calculateNcpArr(
     """Calculates the matrix of NCP from matrix of best fit parameters"""
 
     allNcpForEachMass = []
+    allFseForEachMass = []
     for i in range(len(arrBestFitPars)):
-        ncpForEachMass = calculateNcpRow(
+        ncpForEachMass, fseForEachMass = calculateNcpRow(
             arrBestFitPars[i],
             ySpacesForEachMass[i],
             resolutionPars[i],
@@ -448,10 +450,12 @@ def calculateNcpArr(
         )
 
         allNcpForEachMass.append(ncpForEachMass)
+        allFseForEachMass.append(fseForEachMass)
 
     allNcpForEachMass = np.array(allNcpForEachMass)
     allNcpTotal = np.sum(allNcpForEachMass, axis=1)
-    return allNcpForEachMass, allNcpTotal
+    allFseForEachMass = np.array(allFseForEachMass)
+    return allNcpForEachMass, allNcpTotal, allFseForEachMass
 
 
 def calculateNcpRow(
@@ -461,12 +465,12 @@ def calculateNcpRow(
     output: row shape with the ncpTotal for each mass"""
 
     if np.all(initPars == 0):
-        return np.zeros(ySpacesForEachMass.shape)
+        return np.zeros(ySpacesForEachMass.shape), np.zeros(ySpacesForEachMass.shape)
 
-    ncpForEachMass, ncpTotal = calculateNcpSpec(
+    ncpForEachMass, ncpTotal, fseForEachMass = calculateNcpSpec(
         ic, initPars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays
     )
-    return ncpForEachMass
+    return ncpForEachMass, fseForEachMass
 
 
 def createNcpWorkspaces(ncpForEachMass, ncpTotal, ws, ic):
@@ -508,6 +512,31 @@ def createNcpWorkspaces(ncpForEachMass, ncpTotal, ws, ic):
 
     return wsTotNCPSum, wsMNCPSum
 
+def createFseWorkspaces(fseForEachMass, ws, ic):
+    # Need to rearrage array of yspaces into seperate arrays for each mass
+    fseForEachMass = switchFirstTwoAxis(fseForEachMass)
+
+    # Use ws dataX to match with histogram data
+    dataX = ws.extractX()[
+        :, : fseForEachMass.shape[2]
+    ]  # Make dataX match ncp shape automatically
+    assert (
+        fseForEachMass[0].shape == dataX.shape
+    ), "DataX and DataY in ws need to be the same shape."
+
+    # Individual ncp workspaces
+    wsMFSESum = []
+    for i, fse_m in enumerate(fseForEachMass):
+        ws_m = CloneWorkspace(ws, OutputWorkspace=ws.name() + "_TOF_FSE_" + str(i))
+        passDataIntoWS(dataX, fse_m, np.zeros_like(dataX), ws_m)
+        MaskDetectors(Workspace=ws_m, WorkspaceIndexList=ic.maskedDetectorIdx)
+        ws_m_sum = SumSpectra(
+            InputWorkspace=ws_m, OutputWorkspace=ws_m.name() + "_Sum"
+        )
+        wsMFSESum.append(ws_m_sum)
+
+    return wsMFSESum
+
 
 def createWS(dataX, dataY, dataE, wsName):
     ws = CreateWorkspace(
@@ -710,7 +739,7 @@ def errorFunction(
 ):
     """Error function to be minimized, operates in TOF space"""
 
-    ncpForEachMass, ncpTotal = calculateNcpSpec(
+    ncpForEachMass, ncpTotal, fseForEachMass = calculateNcpSpec(
         ic, pars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays
     )
 
@@ -752,7 +781,7 @@ def calculateNcpSpec(
 
     ncpForEachMass = intensities * (JOfY + FSE) * E0 * E0 ** (-0.92) * masses / deltaQ
     ncpTotal = np.sum(ncpForEachMass, axis=0)
-    return ncpForEachMass, ncpTotal
+    return ncpForEachMass, ncpTotal, FSE
 
 
 def prepareArraysFromPars(ic, initPars):
diff --git a/mvesuvio/vesuvio_analysis/fit_in_yspace.py b/mvesuvio/vesuvio_analysis/fit_in_yspace.py
index c8ea26e1..81f07180 100644
--- a/mvesuvio/vesuvio_analysis/fit_in_yspace.py
+++ b/mvesuvio/vesuvio_analysis/fit_in_yspace.py
@@ -18,6 +18,9 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF):
 
     wsTOFMass0 = subtractAllMassesExceptFirst(IC, wsTOF, ncpForEachMass)
 
+    if yFitIC.subtractFSE:
+        wsTOFMass0 = Minus(wsTOFMass0, wsTOF.name() + "_TOF_FSE_0", OutputWorkspace=wsTOFMass0.name() + "_FSE")
+
     wsJoY, wsJoYAvg = ySpaceReduction(
         wsTOFMass0, IC.masses[0], yFitIC, ncpForEachMass[:, 0, :]
     )