diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LiquidsReflectometryReduction.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LiquidsReflectometryReduction.py
index 1cf07a5b262690f0c77b3417246447a136ba0102..db08a5e528100ac0f8780f28cfc59ec46d76d0df 100644
--- a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LiquidsReflectometryReduction.py
+++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LiquidsReflectometryReduction.py
@@ -2,7 +2,6 @@
 import time
 import math
 import os
-import numpy
 from mantid.api import *
 from mantid.simpleapi import *
 from mantid.kernel import *
@@ -54,7 +53,7 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
                                                 "TOF range to use")
         self.declareProperty("TofRangeFlag", True,
                              doc="If true, the TOF will be cropped according to the TOF range property")
-        self.declareProperty("QMin", 0.05, doc="Mnimum Q-value")
+        self.declareProperty("QMin", 0.05, doc="Minimum Q-value")
         self.declareProperty("QStep", 0.02, doc="Step size in Q. Enter a negative value to get a log scale")
         self.declareProperty("AngleOffset", 0.0, doc="angle offset (degrees)")
         self.declareProperty("AngleOffsetError", 0.0, doc="Angle offset error (degrees)")
@@ -66,6 +65,8 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         self.declareProperty("GeometryCorrectionFlag", False, doc="Use or not the geometry correction")
         self.declareProperty("FrontSlitName", "S1", doc="Name of the front slit")
         self.declareProperty("BackSlitName", "Si", doc="Name of the back slit")
+        self.declareProperty("TOFSteps", 40.0, doc="TOF step size")
+        self.declareProperty("CropFirstAndLastPoints", True, doc="If true, we crop the first and last points")
 
     def PyExec(self):
         # The old reflectivity reduction has an offset between the input
@@ -87,14 +88,7 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         normBackRange = self.getProperty("NormBackgroundPixelRange").value
         normPeakRange = self.getProperty("NormPeakPixelRange").value
 
-        #GENERAL
-        #TODO: Why are there two versions of this?
-        TOFrangeFlag = self.getProperty("TofRangeFlag")
-        if TOFrangeFlag:
-            TOFrange = self.getProperty("TOFRange").value  #microS
-        else:
-            TOFrange = [0, 200000]
-
+        # Get Q range
         qMin = self.getProperty("QMin").value
         qStep = self.getProperty("QStep").value
         if qStep > 0:  #force logarithmic binning
@@ -108,20 +102,19 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         runs = reduce((lambda x, y: '%s+%s' % (x, y)), file_list)
         ws_event_data = Load(Filename=runs)
 
+        # Get the TOF range
+        TOFrangeFlag = self.getProperty("TofRangeFlag")
+        if TOFrangeFlag:
+            TOFrange = self.getProperty("TOFRange").value  #microS
+        else:
+            tof_max = ws_event_data.getTofMax()
+            TOFrange = [0, tof_max]
+
         # Number of pixels in each direction
         #TODO: revisit this when we update the IDF
         self.number_of_pixels_x = int(ws_event_data.getInstrument().getNumberParameter("number-of-x-pixels")[0])
         self.number_of_pixels_y = int(ws_event_data.getInstrument().getNumberParameter("number-of-y-pixels")[0])
 
-        # get the distance moderator-detector and sample-detector
-        #TODO: get rid of this
-        sample = ws_event_data.getInstrument().getSample()
-        source = ws_event_data.getInstrument().getSource()
-        source_sample_distance = sample.getDistance(source)
-        detector = ws_event_data.getDetector(0)
-        sample_detector_distance = detector.getPos().getZ()
-        source_detector_distance = source_sample_distance + sample_detector_distance
-
         # Get scattering angle theta
         theta = self.calculate_scattering_angle(ws_event_data)
 
@@ -155,8 +148,11 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         AnalysisDataService.remove(str(data_cropped))
         AnalysisDataService.remove(str(norm_summed))
 
+        normalized_data = ConvertToPointData(InputWorkspace=normalized_data,
+                                             OutputWorkspace=str(normalized_data))
+
         # Apply scaling factors
-        self.apply_scaling_factor(normalized_data)
+        normalized_data = self.apply_scaling_factor(normalized_data)
 
         q_workspace = SumSpectra(InputWorkspace = normalized_data)
         q_workspace.getAxis(0).setUnit("MomentumTransfer")
@@ -166,8 +162,15 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         if geometryCorrectionFlag:
             logger.error("The geometry correction for the Q conversion has not been implemented.")
 
-        # create workspace
-        # Q range
+        # Get the distance fromthe moderator to the detector
+        sample = ws_event_data.getInstrument().getSample()
+        source = ws_event_data.getInstrument().getSource()
+        source_sample_distance = sample.getDistance(source)
+        detector = ws_event_data.getDetector(0)
+        sample_detector_distance = detector.getPos().getZ()
+        source_detector_distance = source_sample_distance + sample_detector_distance
+
+        # Convert to Q
         # Use the TOF range to pick the maximum Q, and give it a little extra room.
         h = 6.626e-34  # m^2 kg s^-1
         m = 1.675e-27  # kg
@@ -178,19 +181,16 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         for i in range(len(data_x)):
             data_x[i] = constant / data_x[i]
         q_workspace = SortXAxis(InputWorkspace=q_workspace, OutputWorkspace=str(q_workspace))
-            
-        
+
+        # Cook up a name compatible with the UI for backward compatibility
         _time = int(time.time())
         name_output_ws = self.getPropertyValue("OutputWorkspace")
         name_output_ws = name_output_ws + '_#' + str(_time) + 'ts'
 
         q_rebin = Rebin(InputWorkspace=q_workspace, Params=q_range,
                         OutputWorkspace=name_output_ws)
-        q_rebin = ConvertToPointData(InputWorkspace=q_rebin,
-                                     OutputWorkspace=name_output_ws)
 
         # Crop to non-zero values
-        # Get rid of first and last Q points to avoid edge effects
         data_y = q_rebin.readY(0)
         low_q = None
         high_q = None
@@ -201,19 +201,35 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
                 high_q = len(data_y)-1-i
             if low_q is not None and high_q is not None:
                 break
-            
+
+        crop = self.getProperty("CropFirstAndLastPoints").value
         if low_q is not None and high_q is not None:
+            # Get rid of first and last Q points to avoid edge effects
+            if crop:
+                low_q += 1
+                high_q -= 1
             data_x = q_rebin.readX(0)
             q_rebin = CropWorkspace(InputWorkspace=q_rebin,
                                     OutputWorkspace=str(q_rebin),
-                                    XMin=data_x[low_q+1], XMax=data_x[high_q-1])
+                                    XMin=data_x[low_q], XMax=data_x[high_q])
         else:
             logger.error("Data is all zeros. Check your TOF ranges.")
 
-        #TODO: why do we need this cleanup?
-        # Cleanup data turns values < 1e-12 and values where the error is greater than the value
-        # to 0 +- 1
-        #[final_y_axis, final_y_error_axis] = wks_utility.cleanupData1D(final_y_axis, final_error_axis)
+        # Clean up the workspace for backward compatibility
+        data_y = q_rebin.dataY(0)
+        data_e = q_rebin.dataE(0)
+        # Again for backward compatibility, the first and last points of the 
+        # raw output when not cropping was simply set to 0 += 1.
+        if crop is False:
+            data_y[0] = 0
+            data_e[0] = 1
+            data_y[len(data_y)-1] = 0
+            data_e[len(data_y)-1] = 1
+        # Values < 1e-12 and values where the error is greater than the value are replaced by 0+-1
+        for i in range(len(data_y)):
+            if data_y[i] < 1e-12 or data_e[i]>data_y[i]:
+                data_y[i]=0.0
+                data_e[i]=1.0
 
         #TODO: remove this, which we use during development to make sure we don't leave trash
         logger.information(str(AnalysisDataService.getObjectNames()))
@@ -275,7 +291,8 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         x_min = int(low_res_range[0]) + offset
         x_max = int(low_res_range[1]) + offset
 
-        if peak_range[0] > background_range[0] and peak_range[1] < background_range[1]:
+        left_bck = None
+        if peak_min > bck_min:
             left_bck = RefRoi(InputWorkspace=workspace, IntegrateY=False,
                               NXPixel=self.number_of_pixels_x,
                               NYPixel=self.number_of_pixels_y,
@@ -284,8 +301,11 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
                               XPixelMax=x_max,
                               YPixelMin=bck_min,
                               YPixelMax=peak_min - 1,
+                              ErrorWeighting = True,
                               SumPixels=True, NormalizeSum=True)
 
+        right_bck = None
+        if peak_max < bck_max:
             right_bck = RefRoi(InputWorkspace=workspace, IntegrateY=False,
                                NXPixel=self.number_of_pixels_x,
                                NYPixel=self.number_of_pixels_y,
@@ -294,11 +314,15 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
                                XPixelMax=x_max,
                                YPixelMin=peak_max + 1,
                                YPixelMax=bck_max,
+                               ErrorWeighting = True,
                                SumPixels=True, NormalizeSum=True)
+            
+        if right_bck is not None and left_bck is not None:
             average = (left_bck + right_bck) / 2.0
-            # Avoid leaving trash behind
-            AnalysisDataService.remove(str(left_bck))
-            AnalysisDataService.remove(str(right_bck))
+        elif right_bck is not None:
+            average = right_bck
+        elif left_bck is not None:
+            average = left_bck
         else:
             average = RefRoi(InputWorkspace=workspace, IntegrateY=False,
                              NXPixel=self.number_of_pixels_x,
@@ -308,6 +332,7 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
                              XPixelMax=x_max,
                              YPixelMin=bck_min,
                              YPixelMax=bck_max,
+                             ErrorWeighting = True,
                              SumPixels=True, NormalizeSum=True)
         # Integrate over the low-res direction
         workspace = RefRoi(InputWorkspace=workspace, IntegrateY=False,
@@ -318,10 +343,18 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
                            ConvertToQ=False,
                            SumPixels=sum_peak,
                            OutputWorkspace=str(workspace))
+        #TODO Check whether we should multiply by the number of pixels
+        # in the low-res direction
         workspace = Minus(LHSWorkspace=workspace, RHSWorkspace=average,
                           OutputWorkspace=str(workspace))
         # Avoid leaving trash behind
-        AnalysisDataService.remove(str(average))
+        average_name = str(average)
+        if AnalysisDataService.doesExist(str(left_bck)):
+            AnalysisDataService.remove(str(left_bck))
+        if AnalysisDataService.doesExist(str(right_bck)):
+            AnalysisDataService.remove(str(right_bck))
+        if AnalysisDataService.doesExist(average_name):
+            AnalysisDataService.remove(average_name)
         return workspace
 
     def process_data(self, workspace, tof_range, crop_low_res, low_res_range,
@@ -329,9 +362,10 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
         """
             Common processing for both sample data and normalization.
         """
-        # Rebin normalization
-        #TODO: Check why they chose those hard-coded values
-        workspace = Rebin(InputWorkspace=workspace, Params=[0, 40, 200000], 
+        # Rebin TOF axis
+        tof_max = workspace.getTofMax()
+        tof_step = self.getProperty("TOFSteps").value
+        workspace = Rebin(InputWorkspace=workspace, Params=[0, tof_step, tof_max], 
                           PreserveEvents=False, OutputWorkspace="%s_histo" % str(workspace))
 
         # Crop TOF range
@@ -493,8 +527,7 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
             # Should it not be done for the center of the bin?
             for i in range(len(norm_value)):
                 norm_value[i] = norm_tof[i] * b + a
-                #TODO: Verify with REFL team that they really meant the following line
-                norm_error[i] = norm_tof[i] * b_error + a_error
+                norm_error[i] = math.sqrt(a_error*a_error + norm_tof[i] * norm_tof[i] * b_error * b_error)
 
             workspace = Divide(LHSWorkspace=workspace,
                                RHSWorkspace=normalization,
diff --git a/Code/Mantid/Framework/WorkflowAlgorithms/src/RefReduction.cpp b/Code/Mantid/Framework/WorkflowAlgorithms/src/RefReduction.cpp
index c4f1a9b4ccc9ba04a40190f297a759bb3b1ed778..a861fefdef291d4b5270482e8d423fd0e2e30be0 100644
--- a/Code/Mantid/Framework/WorkflowAlgorithms/src/RefReduction.cpp
+++ b/Code/Mantid/Framework/WorkflowAlgorithms/src/RefReduction.cpp
@@ -443,7 +443,7 @@ MatrixWorkspace_sptr RefReduction::processNormalization() {
   refAlg->setProperty("ConvertToQ", false);
   refAlg->setProperty("SumPixels", true);
   refAlg->setProperty("NormalizeSum", true);
-
+  refAlg->setProperty("AverageOverIntegratedAxis", integrateY);
   refAlg->setProperty("YPixelMin", ymin);
   refAlg->setProperty("YPixelMax", ymax);
   refAlg->setProperty("XPixelMin", xmin);
@@ -712,7 +712,7 @@ MatrixWorkspace_sptr RefReduction::subtractBackground(
     leftAlg->setProperty("ConvertToQ", false);
     leftAlg->setProperty("SumPixels", true);
     leftAlg->setProperty("NormalizeSum", true);
-
+    leftAlg->setProperty("AverageOverIntegratedAxis", integrateY);
     leftAlg->setProperty("YPixelMin", ymin);
     leftAlg->setProperty("YPixelMax", ymax);
     leftAlg->setProperty("XPixelMin", xmin);
@@ -737,7 +737,7 @@ MatrixWorkspace_sptr RefReduction::subtractBackground(
     rightAlg->setProperty("ConvertToQ", false);
     rightAlg->setProperty("SumPixels", true);
     rightAlg->setProperty("NormalizeSum", true);
-
+    rightAlg->setProperty("AverageOverIntegratedAxis", integrateY);
     rightAlg->setProperty("YPixelMin", ymin);
     rightAlg->setProperty("YPixelMax", ymax);
     rightAlg->setProperty("XPixelMin", xmin);
@@ -778,7 +778,7 @@ MatrixWorkspace_sptr RefReduction::subtractBackground(
     refAlg->setProperty("ConvertToQ", false);
     refAlg->setProperty("SumPixels", true);
     refAlg->setProperty("NormalizeSum", true);
-
+    refAlg->setProperty("AverageOverIntegratedAxis", integrateY);
     refAlg->setProperty("YPixelMin", ymin);
     refAlg->setProperty("YPixelMax", ymax);
     refAlg->setProperty("XPixelMin", xmin);
diff --git a/Code/Mantid/Framework/WorkflowAlgorithms/src/RefRoi.cpp b/Code/Mantid/Framework/WorkflowAlgorithms/src/RefRoi.cpp
index 50be81a9ad8d7f491122765582715d7b7c88d232..20dbd035fe1cb61257fdaee79f6bb7a826764163 100644
--- a/Code/Mantid/Framework/WorkflowAlgorithms/src/RefRoi.cpp
+++ b/Code/Mantid/Framework/WorkflowAlgorithms/src/RefRoi.cpp
@@ -46,7 +46,14 @@ void RefRoi::init() {
   declareProperty(
       "NormalizeSum", false,
       "If true, and SumPixels is true, the"
-      "resulting histogram will be divided by the number of pixels in the ROI");
+      " resulting histogram will be divided by the number of pixels in the ROI");
+  declareProperty(
+      "AverageOverIntegratedAxis", false,
+      "If true, and SumPixels and NormalizeSum are true, the"
+      " resulting histogram will also be divided by the number of pixels integrated over");
+  declareProperty(
+      "ErrorWeighting", false,
+      "If true, error weighting will be used when normalizing");
   declareProperty(
       "IntegrateY", true,
       "If true, the Y direction will be"
@@ -102,6 +109,7 @@ void RefRoi::extract2D() {
   bool integrate_y = getProperty("IntegrateY");
   bool sum_pixels = getProperty("SumPixels");
   bool normalize = getProperty("NormalizeSum");
+  bool error_weighting = getProperty("ErrorWeighting");
 
   int nHisto = integrate_y ? m_nXPixel : m_nYPixel;
   int xmin = integrate_y ? 0 : m_xMin;
@@ -149,44 +157,74 @@ void RefRoi::extract2D() {
     XOut0 = inputWS->readX(0);
   }
 
-  //PARALLEL_FOR2(outputWS, inputWS)
-  for (int i = xmin; i <= xmax; i++) {
-    //PARALLEL_START_INTERUPT_REGION
-    for (int j = ymin; j <= ymax; j++) {
-      int index = m_nYPixel * i + j;
+  // Make sure the inner loop is always the one we integrate over
+  int main_axis_min = integrate_y ? xmin : ymin;
+  int main_axis_max = integrate_y ? xmax : ymax;
+  int integrated_axis_min = integrate_y ? ymin : xmin;
+  int integrated_axis_max = integrate_y ? ymax : xmax;
+
+  for (int i = main_axis_min; i <= main_axis_max; i++) {
+    size_t output_index = i;
+    if (sum_pixels)
+      output_index = 0;
+
+    MantidVec &YOut = outputWS->dataY(output_index);
+    MantidVec &EOut = outputWS->dataE(output_index);
+    MantidVec signal_vector(YOut.size(), 0.0);
+    MantidVec error_vector(YOut.size(), 0.0);
+
+    for (int j = integrated_axis_min; j <= integrated_axis_max; j++) {
+      int index = integrate_y ? m_nYPixel * i + j : m_nYPixel * j + i;
       const MantidVec &YIn = inputWS->readY(index);
       const MantidVec &EIn = inputWS->readE(index);
 
-      size_t output_index = integrate_y ? i : j;
-      if (sum_pixels)
-        output_index = 0;
-
-      MantidVec &YOut = outputWS->dataY(output_index);
-      MantidVec &EOut = outputWS->dataE(output_index);
+      for (size_t t = 0; t < YOut.size(); t++) {
+        size_t t_index = convert_to_q ? YOut.size() - 1 - t : t;
+        if (sum_pixels && normalize && error_weighting) {
+          signal_vector[t] += YIn[t_index];
+          error_vector[t] += EIn[t_index] * EIn[t_index];
+        } else {
+          YOut[t] += YIn[t_index];
+          EOut[t] += EIn[t_index] * EIn[t_index];
+        }
+      }
+    }
 
+    if (sum_pixels && normalize && error_weighting) {
       for (size_t t = 0; t < YOut.size(); t++) {
         size_t t_index = convert_to_q ? YOut.size() - 1 - t : t;
-        YOut[t] += YIn[t_index];
-        EOut[t] += EIn[t_index] * EIn[t_index];
+        double error_squared = error_vector[t_index] == 0 ? 1 : error_vector[t_index];
+        YOut[t] += signal_vector[t_index] / error_squared;
+        EOut[t] += 1.0 / error_squared;
       }
     }
-    //PARALLEL_END_INTERUPT_REGION
   }
-  //PARALLEL_CHECK_INTERUPT_REGION
 
-  const int n_pixels = (xmax - xmin + 1) * (ymax - ymin + 1);
+  // Check whether we want to divide by the number of pixels along
+  // the axis we integrated over.
+  bool average_integrated = getProperty("AverageOverIntegratedAxis");
+  double n_integrated = 1.0;
+  if (sum_pixels && normalize && average_integrated) {
+    n_integrated = integrated_axis_max - integrated_axis_min + 1;
+  }
 
   for (int i = 0; i < nHisto; i++) {
     outputWS->dataX(i) = XOut0;
     MantidVec &YOut = outputWS->dataY(i);
     MantidVec &EOut = outputWS->dataE(i);
     for (size_t t = 0; t < EOut.size(); t++) {
-      EOut[t] = sqrt(EOut[t]);
-
       if (sum_pixels && normalize) {
-        YOut[t] = YOut[t] / n_pixels;
-        EOut[t] = EOut[t] / n_pixels;
-      }
+        if (error_weighting) {
+          YOut[t] = YOut[t] / EOut[t] / n_integrated;
+          EOut[t] = sqrt(1.0/EOut[t]) / n_integrated;
+        } else {
+          EOut[t] = sqrt(EOut[t]);
+          YOut[t] = YOut[t] / (main_axis_max - main_axis_min + 1) / n_integrated;
+          EOut[t] = EOut[t] / (main_axis_max - main_axis_min + 1) / n_integrated;
+        }
+      } else {
+        EOut[t] = sqrt(EOut[t]);
+      };
     }
   }
 
diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/LiquidsReflectometryReductionTest.py b/Code/Mantid/Testing/SystemTests/tests/analysis/LiquidsReflectometryReductionTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..118f1a63c428ddc95043170ae6779182e4ceb76e
--- /dev/null
+++ b/Code/Mantid/Testing/SystemTests/tests/analysis/LiquidsReflectometryReductionTest.py
@@ -0,0 +1,49 @@
+#pylint: disable=no-init
+import stresstesting
+from mantid import *
+
+from mantid.simpleapi import *
+
+class LiquidsReflectometryReductionTest(stresstesting.MantidStressTest):
+    def runTest(self):
+        #TODO: The reduction algorithm should not require an absolute path
+        scaling_factor_file = FileFinder.getFullPath("directBeamDatabaseFall2014_IPTS_11601_2.cfg")
+        
+        LiquidsReflectometryReduction(RunNumbers=[119814],
+                                      NormalizationRunNumber=119690,
+                                      SignalPeakPixelRange=[154, 166],
+                                      SubtractSignalBackground=True,
+                                      SignalBackgroundPixelRange=[151, 169],
+                                      NormFlag=True,
+                                      NormPeakPixelRange=[154, 160],
+                                      NormBackgroundPixelRange=[151, 163],
+                                      SubtractNormBackground=True,
+                                      LowResDataAxisPixelRangeFlag=True,
+                                      LowResDataAxisPixelRange=[99, 158],
+                                      LowResNormAxisPixelRangeFlag=True,
+                                      LowResNormAxisPixelRange=[98, 158],
+                                      TOFRange=[29623.0, 42438.0],
+                                      IncidentMediumSelected='2InDiamSi',
+                                      GeometryCorrectionFlag=False,
+                                      QMin=0.005,
+                                      QStep=0.01,
+                                      AngleOffset=0.009,
+                                      AngleOffsetError=0.001,
+                                      ScalingFactorFile=scaling_factor_file,
+                                      SlitsWidthFlag=True,
+                                      CropFirstAndLastPoints=False,
+                                      OutputWorkspace='reflectivity_119814')
+
+    def validate(self):
+        # Be more tolerant with the output.
+        self.tolerance = 0.01
+
+        # Skip the first point so we don't have to have a big tolerance
+        data_y = mtd["reflectivity_119814"].dataY(0)
+        data_y[1] = 0.631281639115562
+        
+        self.disableChecking.append('Instrument')
+        self.disableChecking.append('Sample')
+        self.disableChecking.append('SpectraMap')
+        self.disableChecking.append('Axes')
+        return "reflectivity_119814", 'REFL_119814_combined_data.nxs'
diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/LiquidsReflectometryReductionWithBackgroundTest.py b/Code/Mantid/Testing/SystemTests/tests/analysis/LiquidsReflectometryReductionWithBackgroundTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d09d62bec2c46703814a59c927399e289099fa5
--- /dev/null
+++ b/Code/Mantid/Testing/SystemTests/tests/analysis/LiquidsReflectometryReductionWithBackgroundTest.py
@@ -0,0 +1,49 @@
+#pylint: disable=no-init
+import stresstesting
+from mantid import *
+from mantid.simpleapi import *
+
+class LiquidsReflectometryReductionWithBackgroundTest(stresstesting.MantidStressTest):
+    def runTest(self):
+        #TODO: The reduction algorithm should not require an absolute path
+        scaling_factor_file = FileFinder.getFullPath("directBeamDatabaseFall2014_IPTS_11601_2.cfg")
+        
+        LiquidsReflectometryReduction(RunNumbers=[119816],
+                                      NormalizationRunNumber=119692,
+                                      SignalPeakPixelRange=[155, 165],
+                                      SubtractSignalBackground=True,
+                                      SignalBackgroundPixelRange=[146, 165],
+                                      NormFlag=True,
+                                      NormPeakPixelRange=[154, 162],
+                                      NormBackgroundPixelRange=[151, 165],
+                                      SubtractNormBackground=True,
+                                      LowResDataAxisPixelRangeFlag=True,
+                                      LowResDataAxisPixelRange=[99, 158],
+                                      LowResNormAxisPixelRangeFlag=True,
+                                      LowResNormAxisPixelRange=[118, 137],
+                                      TOFRange=[9610, 22425],
+                                      IncidentMediumSelected='2InDiamSi',
+                                      GeometryCorrectionFlag=False,
+                                      QMin=0.005,
+                                      QStep=0.01,
+                                      AngleOffset=0.009,
+                                      AngleOffsetError=0.001,
+                                      ScalingFactorFile=scaling_factor_file,
+                                      SlitsWidthFlag=True,
+                                      CropFirstAndLastPoints=False,
+                                      OutputWorkspace='reflectivity_119816')
+
+    def validate(self):
+        # Be more tolerant with the output.
+        self.tolerance = 0.0002
+
+        # Skip the first point so we don't have to have a big tolerance
+        data_y = mtd["reflectivity_119816"].dataY(0)
+        data_y[1] = 0.00499601750282373
+        
+        self.disableChecking.append('Instrument')
+        self.disableChecking.append('Sample')
+        self.disableChecking.append('SpectraMap')
+        self.disableChecking.append('Axes')
+        return "reflectivity_119816", 'REFL_119816.nxs'
+
diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/REFLReduction.py b/Code/Mantid/Testing/SystemTests/tests/analysis/REFLReduction.py
index 85669b4377e89abe1939000eb7ff547438962d4a..57c9b9c7b02505aaa022e88ed3c9dfd830b9429d 100644
--- a/Code/Mantid/Testing/SystemTests/tests/analysis/REFLReduction.py
+++ b/Code/Mantid/Testing/SystemTests/tests/analysis/REFLReduction.py
@@ -1,7 +1,6 @@
 #pylint: disable=no-init
 import stresstesting
 from mantid import *
-
 from mantid.simpleapi import *
 
 class REFLReduction(stresstesting.MantidStressTest):
@@ -9,8 +8,7 @@ class REFLReduction(stresstesting.MantidStressTest):
         #TODO: The reduction algorithm should not require an absolute path
         scaling_factor_file = FileFinder.getFullPath("directBeamDatabaseFall2014_IPTS_11601_2.cfg")
         
-        #RefLReduction(RunNumbers=[119814],
-        LiquidsReflectometryReduction(RunNumbers=[119814],
+        RefLReduction(RunNumbers=[119814],
                       NormalizationRunNumber=119690,
                       SignalPeakPixelRange=[154, 166],
                       SubtractSignalBackground=True,
@@ -36,11 +34,10 @@ class REFLReduction(stresstesting.MantidStressTest):
 
     def validate(self):
         # Be more tolerant with the output.
-        self.tolerance = 0.008
+        self.tolerance = 0.0001
 
         self.disableChecking.append('Instrument')
         self.disableChecking.append('Sample')
         self.disableChecking.append('SpectraMap')
         self.disableChecking.append('Axes')
         return "reflectivity_119814", 'REFL_119814_combined_data.nxs'
-
diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/REFLWithBackground.py b/Code/Mantid/Testing/SystemTests/tests/analysis/REFLWithBackground.py
new file mode 100644
index 0000000000000000000000000000000000000000..17443e9c64b2ad225c6b43ca1906d84ffbe7eb42
--- /dev/null
+++ b/Code/Mantid/Testing/SystemTests/tests/analysis/REFLWithBackground.py
@@ -0,0 +1,44 @@
+#pylint: disable=no-init
+import stresstesting
+from mantid import *
+from mantid.simpleapi import *
+
+class REFLWithBackground(stresstesting.MantidStressTest):
+    def runTest(self):
+        #TODO: The reduction algorithm should not require an absolute path
+        scaling_factor_file = FileFinder.getFullPath("directBeamDatabaseFall2014_IPTS_11601_2.cfg")
+        
+        RefLReduction(RunNumbers=[119816],
+                      NormalizationRunNumber=119692,
+                      SignalPeakPixelRange=[155, 165],
+                      SubtractSignalBackground=True,
+                      SignalBackgroundPixelRange=[146, 165],
+                      NormFlag=True,
+                      NormPeakPixelRange=[154, 162],
+                      NormBackgroundPixelRange=[151, 165],
+                      SubtractNormBackground=True,
+                      LowResDataAxisPixelRangeFlag=True,
+                      LowResDataAxisPixelRange=[99, 158],
+                      LowResNormAxisPixelRangeFlag=True,
+                      LowResNormAxisPixelRange=[118, 137],
+                      TOFRange=[9610, 22425],
+                      IncidentMediumSelected='2InDiamSi',
+                      GeometryCorrectionFlag=False,
+                      QMin=0.005,
+                      QStep=0.01,
+                      AngleOffset=0.009,
+                      AngleOffsetError=0.001,
+                      ScalingFactorFile=scaling_factor_file,
+                      SlitsWidthFlag=True,
+                      OutputWorkspace='reflectivity_119816')
+
+    def validate(self):
+        # Be more tolerant with the output.
+        self.tolerance = 0.0001
+
+        self.disableChecking.append('Instrument')
+        self.disableChecking.append('Sample')
+        self.disableChecking.append('SpectraMap')
+        self.disableChecking.append('Axes')
+        return "reflectivity_119816", 'REFL_119816.nxs'
+
diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/RefRoi.py b/Code/Mantid/Testing/SystemTests/tests/analysis/RefRoi.py
new file mode 100644
index 0000000000000000000000000000000000000000..a8aee4424f5e8bac1694d01217740e048f523edc
--- /dev/null
+++ b/Code/Mantid/Testing/SystemTests/tests/analysis/RefRoi.py
@@ -0,0 +1,22 @@
+#pylint: disable=no-init
+import stresstesting
+from mantid import *
+from mantid.simpleapi import *
+
+class RefRoiTest(stresstesting.MantidStressTest):
+    def runTest(self):
+        ws=Load(Filename="REF_L_119814")
+        ws = Integration(InputWorkspace=ws)
+        roi = RefRoi(InputWorkspace=ws,
+                     NXPixel=256, NYPixel=304,
+                     IntegrateY=False, ConvertToQ=False)
+        roi = Transpose(InputWorkspace=roi)
+
+    def validate(self):
+        self.disableChecking.append('Instrument')
+        self.disableChecking.append('Sample')
+        self.disableChecking.append('SpectraMap')
+        self.disableChecking.append('Axes')
+        return "roi", 'REFL_119814_roi_peak.nxs'
+
+
diff --git a/Code/Mantid/scripts/reduction/instruments/reflectometer/wks_utility.py b/Code/Mantid/scripts/reduction/instruments/reflectometer/wks_utility.py
index 37b9f0d1085a257f4d648f0d0742272909aa3c86..707396dfa9f5ecb92dd68ae2531a09c04ae7cea4 100644
--- a/Code/Mantid/scripts/reduction/instruments/reflectometer/wks_utility.py
+++ b/Code/Mantid/scripts/reduction/instruments/reflectometer/wks_utility.py
@@ -1167,8 +1167,6 @@ def substractBackground(tof_axis, y_axis, y_error_axis,
     return [final_y_axis, final_y_error_axis]
 
 def weightedMean(data_array, error_array, error_0):
-    import numpy
-    return [numpy.average(data_array)/len(data_array), numpy.average(error_array)/len(error_array)]
 
     sz = len(data_array)