diff --git a/Testing/Data/SystemTest/MUT00053591.NXS.md5 b/Testing/Data/SystemTest/MUT00053591.NXS.md5
new file mode 100644
index 0000000000000000000000000000000000000000..778206487c758a2ac3226c290c3068506ae08a63
--- /dev/null
+++ b/Testing/Data/SystemTest/MUT00053591.NXS.md5
@@ -0,0 +1 @@
+125543727e1387609d3bdc22a78890a5
\ No newline at end of file
diff --git a/Testing/SystemTests/tests/analysis/MuonKerenFittingTest.py b/Testing/SystemTests/tests/analysis/MuonKerenFittingTest.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..39fc5d6f87be5c8c32b14e628a4f11edf19e44f0 100644
--- a/Testing/SystemTests/tests/analysis/MuonKerenFittingTest.py
+++ b/Testing/SystemTests/tests/analysis/MuonKerenFittingTest.py
@@ -0,0 +1,39 @@
+#pylint: disable=no-init,attribute-defined-outside-init
+import stresstesting
+from mantid.simpleapi import *
+
+class MuonKerenFittingTest(stresstesting.MantidStressTest):
+    '''Tests the Keren fitting function on a real workspace, to check results vs. WiMDA'''
+
+    def runTest(self):
+        # Load dataset
+        MUT53591 = LoadMuonNexus(Filename='MUT00053591.nxs', DetectorGroupingTable='gp')
+
+	# Process like MuonAnalysis interface would
+	processed = MuonProcess(InputWorkspace='MUT53591', Mode='Combined', SummedPeriodSet='1',\
+	ApplyDeadTimeCorrection=False, DetectorGroupingTable='gp', LoadedTimeZero=0, TimeZero=0,\
+	Xmin=0.08, Xmax=10.0, OutputType="PairAsymmetry", PairFirstIndex="0", PairSecondIndex="1",\
+	Alpha=1.0)
+
+	# Fit the Keren function to the data
+	func = "name=FlatBackground,A0=0.1;name=Keren,A=0.1,Delta=0.2,Field=18,Fluct=0.2"
+	Fit(InputWorkspace='processed', Function=func, Output='out', CreateOutput=True)
+
+	# Get fitted parameters
+	params = mtd['out_Parameters']
+	A0 = params.cell(0,1)
+	A = params.cell(1,1)
+	Delta = params.cell(2,1)
+	Field = params.cell(3,1)
+	Fluct = params.cell(4,1)
+	Chisq = params.cell(5,1)
+
+	# Check that params are within the errors of those obtained in WiMDA
+	self.assertTrue(Chisq < 1.1, "Fitted chi-square too large") 
+	self.assertDelta(A0, 0.1623, 0.0046, "Fitted A0 outside errors")
+	self.assertDelta(A, 0.0389, 0.0040, "Fitted A outside errors")
+	self.assertDelta(Delta, 0.96, 0.11, "Fitted Delta outside errors")
+	self.assertDelta(Field, 20.0, 1.0, "Fitted Field outside errors")
+	self.assertDelta(Fluct, 0.1, 0.01, "Fitted Fluct outside errors")
+
+