diff --git a/Framework/Algorithms/src/CopySample.cpp b/Framework/Algorithms/src/CopySample.cpp
index 6bbc07bd5455c981b908aa80a3079b641de90419..2b28594aa3d928ee2f51a0997e6ac57cbc7ac59e 100644
--- a/Framework/Algorithms/src/CopySample.cpp
+++ b/Framework/Algorithms/src/CopySample.cpp
@@ -80,8 +80,8 @@ void CopySample::exec() {
 
   Sample sample;
   // get input sample
-  IMDEventWorkspace_const_sptr inMDWS =
-      boost::dynamic_pointer_cast<const IMDEventWorkspace>(inWS);
+  MultipleExperimentInfos_sptr inMDWS =
+      boost::dynamic_pointer_cast<MultipleExperimentInfos>(inWS);
   if (inMDWS != nullptr) // it is an MD workspace
   {
     int inputSampleNumber = getProperty("MDInputSampleNumber");
@@ -116,9 +116,8 @@ void CopySample::exec() {
   bool copyOrientation = getProperty("CopyOrientationOnly");
 
   // Sample copy;
-
-  IMDEventWorkspace_sptr outMDWS =
-      boost::dynamic_pointer_cast<IMDEventWorkspace>(outWS);
+  MultipleExperimentInfos_sptr outMDWS =
+      boost::dynamic_pointer_cast<MultipleExperimentInfos>(outWS);
   if (outMDWS != nullptr) {
     int outputSampleNumber = getProperty("MDOutputSampleNumber");
     if ((outputSampleNumber == EMPTY_INT()) ||
diff --git a/Framework/Crystal/src/SetUB.cpp b/Framework/Crystal/src/SetUB.cpp
index c2f3a5c13cbb8061c5e980d5b4f07a0c143e2c28..bf8baf565b6c42d9ee4417bb37311c5845a4c90f 100644
--- a/Framework/Crystal/src/SetUB.cpp
+++ b/Framework/Crystal/src/SetUB.cpp
@@ -110,8 +110,8 @@ void SetUB::exec() {
   Workspace_sptr ws = this->getProperty("Workspace");
 
   // Sample copy;
-  IMDEventWorkspace_sptr mdws =
-      boost::dynamic_pointer_cast<IMDEventWorkspace>(ws);
+  MultipleExperimentInfos_sptr mdws =
+      boost::dynamic_pointer_cast<MultipleExperimentInfos>(ws);
   if (mdws != nullptr) {
     int sampleNumber = getProperty("MDSampleNumber");
     if ((sampleNumber == EMPTY_INT()) ||
diff --git a/Framework/PythonInterface/plugins/algorithms/ConvertWANDSCDtoQ.py b/Framework/PythonInterface/plugins/algorithms/ConvertWANDSCDtoQ.py
new file mode 100644
index 0000000000000000000000000000000000000000..42a5df51e70a544c00209d86462da58095c6c957
--- /dev/null
+++ b/Framework/PythonInterface/plugins/algorithms/ConvertWANDSCDtoQ.py
@@ -0,0 +1,327 @@
+from __future__ import (absolute_import, division, print_function)
+from mantid.api import (PythonAlgorithm, AlgorithmFactory,
+                        PropertyMode, WorkspaceProperty, Progress,
+                        IMDHistoWorkspaceProperty, mtd)
+from mantid.kernel import Direction, FloatArrayProperty, FloatArrayLengthValidator, StringListValidator, FloatBoundedValidator
+from mantid import logger
+import numpy as np
+
+
+class ConvertWANDSCDtoQ(PythonAlgorithm):
+
+    def category(self):
+        return 'DataHandling\\Nexus'
+
+    def seeAlso(self):
+        return [ "LoadWANDSCD" ]
+
+    def name(self):
+        return 'ConvertWANDSCDtoQ'
+
+    def summary(self):
+        return 'Convert the output of LoadWANDSCD to Q or HKL'
+
+    def PyInit(self):
+        self.declareProperty(IMDHistoWorkspaceProperty("InputWorkspace", "",
+                                                       optional=PropertyMode.Mandatory,
+                                                       direction=Direction.Input),
+                             "Input Workspace")
+        self.declareProperty(IMDHistoWorkspaceProperty("NormalisationWorkspace", "",
+                                                       optional=PropertyMode.Optional,
+                                                       direction=Direction.Input),
+                             "Workspace to use for normalisation")
+        self.declareProperty(WorkspaceProperty("UBWorkspace", "",
+                                               optional=PropertyMode.Optional,
+                                               direction=Direction.Input),
+                             "Workspace containing the UB matrix to use")
+        self.declareProperty("Wavelength", 1.488, validator=FloatBoundedValidator(0.0), doc="Wavelength to set the workspace")
+        self.declareProperty('NormaliseBy', 'Monitor', StringListValidator(['None', 'Time', 'Monitor']),
+                             "Normalise to monitor, time or None.")
+        self.declareProperty('Frame', 'Q_sample', StringListValidator(['Q_sample', 'HKL']),
+                             "Selects Q-dimensions of the output workspace")
+        self.declareProperty(FloatArrayProperty("Uproj", [1,0,0], FloatArrayLengthValidator(3), direction=Direction.Input),
+                             "Defines the first projection vector of the target Q coordinate system in HKL mode")
+        self.declareProperty(FloatArrayProperty("Vproj", [0,1,0], FloatArrayLengthValidator(3), direction=Direction.Input),
+                             "Defines the second projection vector of the target Q coordinate system in HKL mode")
+        self.declareProperty(FloatArrayProperty("Wproj", [0,0,1], FloatArrayLengthValidator(3), direction=Direction.Input),
+                             "Defines the third projection vector of the target Q coordinate system in HKL mode")
+        self.declareProperty(FloatArrayProperty("BinningDim0", [-8.02,8.02,401], FloatArrayLengthValidator(3), direction=Direction.Input),
+                             "Binning parameters for the 0th dimension. Enter it as a"
+                             "comma-separated list of values with the"
+                             "format: 'minimum,maximum,number_of_bins'.")
+        self.declareProperty(FloatArrayProperty("BinningDim1", [-0.82,0.82,41], FloatArrayLengthValidator(3), direction=Direction.Input),
+                             "Binning parameters for the 1st dimension. Enter it as a"
+                             "comma-separated list of values with the"
+                             "format: 'minimum,maximum,number_of_bins'.")
+        self.declareProperty(FloatArrayProperty("BinningDim2", [-8.02,8.02,401], FloatArrayLengthValidator(3), direction=Direction.Input),
+                             "Binning parameters for the 2nd dimension. Enter it as a"
+                             "comma-separated list of values with the"
+                             "format: 'minimum,maximum,number_of_bins'.")
+        self.declareProperty('KeepTemporaryWorkspaces', False,
+                             "If True the normalization and data workspaces in addition to the normalized data will be outputted")
+        self.declareProperty(WorkspaceProperty("OutputWorkspace", "",
+                                               optional=PropertyMode.Mandatory,
+                                               direction=Direction.Output),
+                             "Output Workspace")
+
+    def validateInputs(self):
+        issues = dict()
+
+        inWS = self.getProperty("InputWorkspace").value
+
+        if inWS.getNumDims() != 3:
+            issues["InputWorkspace"] = "InputWorkspace has wrong number of dimensions, need 3"
+            return issues
+
+        d0 = inWS.getDimension(0)
+        d1 = inWS.getDimension(1)
+        d2 = inWS.getDimension(2)
+        number_of_runs = d2.getNBins()
+
+        if (d0.name is not 'y' or d1.name is not 'x' or d2.name != 'scanIndex'):
+            issues["InputWorkspace"] = "InputWorkspace has wrong dimensions"
+            return issues
+
+        if inWS.getNumExperimentInfo() == 0:
+            issues["InputWorkspace"] = "InputWorkspace is missing ExperimentInfo"
+            return issues
+
+        # Check that all logs are there and are of correct length
+        run = inWS.getExperimentInfo(0).run()
+        for prop in ['duration', 'monitor_count', 's1']:
+            if run.hasProperty(prop):
+                p = run.getProperty(prop).value
+                if np.size(p) != number_of_runs:
+                    issues["InputWorkspace"] = "log {} is of incorrect length".format(prop)
+            else:
+                issues["InputWorkspace"] = "missing log {}".format(prop)
+
+        for prop in ['azimuthal', 'twotheta']:
+            if run.hasProperty(prop):
+                p = run.getProperty(prop).value
+                if np.size(p) != d0.getNBins()*d1.getNBins():
+                    issues["InputWorkspace"] = "log {} is of incorrect length".format(prop)
+            else:
+                issues["InputWorkspace"] = "missing log {}".format(prop)
+
+        normWS = self.getProperty("NormalisationWorkspace").value
+
+        if normWS:
+            nd0 = normWS.getDimension(0)
+            nd1 = normWS.getDimension(1)
+            nd2 = normWS.getDimension(2)
+            if (nd0.name != d0.name or nd0.getNBins() != d0.getNBins()
+               or nd1.name != d1.name or nd1.getNBins() != d1.getNBins()
+               or nd2.name != d2.name):
+                issues["NormalisationWorkspace"] = "NormalisationWorkspace dimensions are not compatible with InputWorkspace"
+
+        ubWS = self.getProperty("UBWorkspace").value
+        if ubWS:
+            try:
+                sample = ubWS.sample()
+            except AttributeError:
+                sample = ubWS.getExperimentInfo(0).sample()
+            if not sample.hasOrientedLattice():
+                issues["UBWorkspace"] = "UBWorkspace does not has an OrientedLattice"
+        else:
+            if self.getProperty("Frame").value == 'HKL':
+                if not inWS.getExperimentInfo(0).sample().hasOrientedLattice():
+                    issues["Frame"] = "HKL selected but neither an UBWorkspace workspace was provided or " \
+                                      "the InputWorkspace has an OrientedLattice"
+
+        return issues
+
+    def PyExec(self):
+        inWS = self.getProperty("InputWorkspace").value
+        normWS = self.getProperty("NormalisationWorkspace").value
+        _norm = bool(normWS)
+
+        dim0_min, dim0_max, dim0_bins = self.getProperty('BinningDim0').value
+        dim1_min, dim1_max, dim1_bins = self.getProperty('BinningDim1').value
+        dim2_min, dim2_max, dim2_bins = self.getProperty('BinningDim2').value
+        dim0_bins = int(dim0_bins)
+        dim1_bins = int(dim1_bins)
+        dim2_bins = int(dim2_bins)
+        dim0_bin_size = (dim0_max-dim0_min)/dim0_bins
+        dim1_bin_size = (dim1_max-dim1_min)/dim1_bins
+        dim2_bin_size = (dim2_max-dim2_min)/dim2_bins
+
+        data_array = inWS.getSignalArray() # getSignalArray returns a F_CONTIGUOUS view of the signal array
+
+        number_of_runs = data_array.shape[2]
+
+        progress = Progress(self, 0.0, 1.0, number_of_runs+4)
+
+        # Get rotation array
+        s1 = np.deg2rad(inWS.getExperimentInfo(0).run().getProperty('s1').value)
+
+        normaliseBy = self.getProperty("NormaliseBy").value
+        if normaliseBy == "Monitor":
+            scale = np.asarray(inWS.getExperimentInfo(0).run().getProperty('monitor_count').value)
+        elif normaliseBy == "Time":
+            scale = np.asarray(inWS.getExperimentInfo(0).run().getProperty('duration').value)
+        else:
+            scale = np.ones(number_of_runs)
+
+        if _norm:
+            if normaliseBy == "Monitor":
+                norm_scale = np.sum(normWS.getExperimentInfo(0).run().getProperty('monitor_count').value)
+            elif normaliseBy == "Time":
+                norm_scale = np.sum(normWS.getExperimentInfo(0).run().getProperty('duration').value)
+            else:
+                norm_scale = 1.
+            norm_array = normWS.getSignalArray().sum(axis=2)
+
+        W = np.eye(3)
+        UBW = np.eye(3)
+        if self.getProperty("Frame").value == 'HKL':
+            W[:,0] = self.getProperty('Uproj').value
+            W[:,1] = self.getProperty('Vproj').value
+            W[:,2] = self.getProperty('Wproj').value
+            ubWS = self.getProperty("UBWorkspace").value
+            if ubWS:
+                try:
+                    ol = ubWS.sample().getOrientedLattice()
+                except AttributeError:
+                    ol = ubWS.getExperimentInfo(0).sample().getOrientedLattice()
+                logger.notice("Using UB matrix from {} with {}".format(ubWS.name(), ol))
+            else:
+                ol = inWS.getExperimentInfo(0).sample().getOrientedLattice()
+                logger.notice("Using UB matrix from {} with {}".format(inWS.name(), ol))
+            UB = ol.getUB()
+            UBW = np.dot(UB, W)
+            char_dict = {0:'0', 1:'{1}', -1:'-{1}'}
+            chars=['H','K','L']
+            names = ['['+','.join(char_dict.get(j, '{0}{1}')
+                                  .format(j,chars[np.argmax(np.abs(W[:,i]))]) for j in W[:,i])+']' for i in range(3)]
+            units = 'in {:.3f} A^-1,in {:.3f} A^-1,in {:.3f} A^-1'.format(ol.qFromHKL(W[0]).norm(),
+                                                                          ol.qFromHKL(W[1]).norm(),
+                                                                          ol.qFromHKL(W[2]).norm())
+            frames = 'HKL,HKL,HKL'
+            k = 1/self.getProperty("Wavelength").value # Not 2pi/wavelength to save dividing by 2pi later
+        else:
+            names = 'Q_sample_x,Q_sample_y,Q_sample_z'
+            units = 'Angstrom^-1,Angstrom^-1,Angstrom^-1'
+            frames = 'QSample,QSample,QSample'
+            k = 2*np.pi/self.getProperty("Wavelength").value
+
+        progress.report('Calculating Qlab for each pixel')
+        polar = np.array(inWS.getExperimentInfo(0).run().getProperty('twotheta').value)
+        azim = np.array(inWS.getExperimentInfo(0).run().getProperty('azimuthal').value)
+        qlab = np.vstack((np.sin(polar)*np.cos(azim),
+                          np.sin(polar)*np.sin(azim),
+                          np.cos(polar) - 1)).T * -k # Kf - Ki(0,0,1)
+
+        progress.report('Calculating Q volume')
+
+        output = np.zeros((dim0_bins+2, dim1_bins+2, dim2_bins+2))
+        outputr = output.ravel()
+
+        output_scale = np.zeros_like(output)
+        output_scaler = output_scale.ravel()
+
+        if _norm:
+            output_norm = np.zeros_like(output)
+            output_normr = output_norm.ravel()
+            output_norm_scale = np.zeros_like(output)
+            output_norm_scaler = output_norm_scale.ravel()
+
+        bin_size = np.array([[dim0_bin_size],
+                             [dim1_bin_size],
+                             [dim2_bin_size]])
+
+        offset = np.array([[dim0_min/dim0_bin_size],
+                           [dim1_min/dim1_bin_size],
+                           [dim2_min/dim2_bin_size]])-0.5
+
+        assert not data_array[:,:,0].flags.owndata
+        assert not data_array[:,:,0].ravel('F').flags.owndata
+        assert data_array[:,:,0].flags.fnc
+
+        for n in range(number_of_runs):
+            R = np.array([[ np.cos(s1[n]), 0, np.sin(s1[n])],
+                          [             0, 1,             0],
+                          [-np.sin(s1[n]), 0, np.cos(s1[n])]])
+            RUBW = np.dot(R,UBW)
+            q = np.round(np.dot(np.linalg.inv(RUBW),qlab.T)/bin_size-offset).astype(np.int)
+            q_index = np.ravel_multi_index(q, (dim0_bins+2, dim1_bins+2, dim2_bins+2), mode='clip')
+            q_uniq, inverse = np.unique(q_index, return_inverse=True)
+            outputr[q_uniq] += np.bincount(inverse, data_array[:,:,n].ravel('F'))
+            output_scaler[q_uniq] += np.bincount(inverse)*scale[n]
+            if _norm:
+                output_normr[q_uniq] += np.bincount(inverse, norm_array.ravel('F'))
+                output_norm_scaler[q_uniq] += np.bincount(inverse)
+
+            progress.report()
+
+        if _norm:
+            output *= output_norm_scale*norm_scale
+            output_norm *= output_scale
+        else:
+            output_norm = output_scale
+
+        if self.getProperty('KeepTemporaryWorkspaces').value:
+            # Create data workspace
+            progress.report('Creating data MDHistoWorkspace')
+            createWS_alg = self.createChildAlgorithm("CreateMDHistoWorkspace", enableLogging=False)
+            createWS_alg.setProperty("SignalInput", output[1:-1,1:-1,1:-1].ravel('F'))
+            createWS_alg.setProperty("ErrorInput", np.sqrt(output[1:-1,1:-1,1:-1].ravel('F')))
+            createWS_alg.setProperty("Dimensionality", 3)
+            createWS_alg.setProperty("Extents", '{},{},{},{},{},{}'.format(dim0_min,dim0_max,dim1_min,dim1_max,dim2_min,dim2_max))
+            createWS_alg.setProperty("NumberOfBins", '{},{},{}'.format(dim0_bins,dim1_bins,dim2_bins))
+            createWS_alg.setProperty("Names", names)
+            createWS_alg.setProperty("Units", units)
+            createWS_alg.setProperty("Frames", frames)
+            createWS_alg.execute()
+            outWS_data = createWS_alg.getProperty("OutputWorkspace").value
+            mtd.addOrReplace(self.getPropertyValue("OutputWorkspace")+'_data', outWS_data)
+
+            # Create normalisation workspace
+            progress.report('Creating norm MDHistoWorkspace')
+            createWS_alg = self.createChildAlgorithm("CreateMDHistoWorkspace", enableLogging=False)
+            createWS_alg.setProperty("SignalInput", output_norm[1:-1,1:-1,1:-1].ravel('F'))
+            createWS_alg.setProperty("ErrorInput", np.sqrt(output_norm[1:-1,1:-1,1:-1].ravel('F')))
+            createWS_alg.setProperty("Dimensionality", 3)
+            createWS_alg.setProperty("Extents", '{},{},{},{},{},{}'.format(dim0_min,dim0_max,dim1_min,dim1_max,dim2_min,dim2_max))
+            createWS_alg.setProperty("NumberOfBins", '{},{},{}'.format(dim0_bins,dim1_bins,dim2_bins))
+            createWS_alg.setProperty("Names", names)
+            createWS_alg.setProperty("Units", units)
+            createWS_alg.setProperty("Frames", frames)
+            createWS_alg.execute()
+            mtd.addOrReplace(self.getPropertyValue("OutputWorkspace")+'_normalization', createWS_alg.getProperty("OutputWorkspace").value)
+
+        output /= output_norm
+
+        progress.report('Creating MDHistoWorkspace')
+        createWS_alg = self.createChildAlgorithm("CreateMDHistoWorkspace", enableLogging=False)
+        createWS_alg.setProperty("SignalInput", output[1:-1,1:-1,1:-1].ravel('F'))
+        createWS_alg.setProperty("ErrorInput", np.sqrt(output[1:-1,1:-1,1:-1].ravel('F')))
+        createWS_alg.setProperty("Dimensionality", 3)
+        createWS_alg.setProperty("Extents", '{},{},{},{},{},{}'.format(dim0_min,dim0_max,dim1_min,dim1_max,dim2_min,dim2_max))
+        createWS_alg.setProperty("NumberOfBins", '{},{},{}'.format(dim0_bins,dim1_bins,dim2_bins))
+        createWS_alg.setProperty("Names", names)
+        createWS_alg.setProperty("Units", units)
+        createWS_alg.setProperty("Frames", frames)
+        createWS_alg.execute()
+        outWS = createWS_alg.getProperty("OutputWorkspace").value
+
+        # Copy experiment infos
+        if inWS.getNumExperimentInfo() > 0:
+            outWS.copyExperimentInfos(inWS)
+
+        outWS.getExperimentInfo(0).run().addProperty('RUBW_MATRIX', list(UBW.flatten()), True)
+        outWS.getExperimentInfo(0).run().addProperty('W_MATRIX', list(W.flatten()), True)
+        try:
+            if outWS.getExperimentInfo(0).sample().hasOrientedLattice():
+                outWS.getExperimentInfo(0).sample().getOrientedLattice().setUB(UB)
+        except NameError:
+            pass
+
+        if self.getProperty('KeepTemporaryWorkspaces').value:
+            outWS_data.copyExperimentInfos(outWS)
+
+        progress.report()
+        self.setProperty("OutputWorkspace", outWS)
+
+
+AlgorithmFactory.subscribe(ConvertWANDSCDtoQ)
diff --git a/Framework/PythonInterface/plugins/algorithms/LoadWANDSCD.py b/Framework/PythonInterface/plugins/algorithms/LoadWANDSCD.py
new file mode 100644
index 0000000000000000000000000000000000000000..13695e77525d3a3dd7dfa756c00405ce79c1a360
--- /dev/null
+++ b/Framework/PythonInterface/plugins/algorithms/LoadWANDSCD.py
@@ -0,0 +1,172 @@
+from __future__ import (absolute_import, division, print_function)
+from mantid.api import PythonAlgorithm, AlgorithmFactory, PropertyMode, WorkspaceProperty, Progress, MultipleFileProperty, FileAction, mtd
+from mantid.kernel import Direction, Property, IntArrayProperty, StringListValidator
+from mantid.simpleapi import LoadEventNexus, RemoveLogs, DeleteWorkspace, ConvertToMD, Rebin, CreateGroupingWorkspace, GroupDetectors, SetUB
+import numpy as np
+import h5py
+import re
+
+
+class LoadWANDSCD(PythonAlgorithm):
+
+    def category(self):
+        return 'DataHandling\\Nexus'
+
+    def seeAlso(self):
+        return [ "ConvertWANDSCDtoQ" ]
+
+    def name(self):
+        return 'LoadWANDSCD'
+
+    def summary(self):
+        return 'Load WAND single crystal data into a detector space vs rotation MDHisto'
+
+    def PyInit(self):
+        self.declareProperty(MultipleFileProperty(name="Filename", action=FileAction.OptionalLoad,
+                                                  extensions=[".nxs.h5"]),
+                             "Files to load")
+        self.declareProperty('IPTS', Property.EMPTY_INT, "IPTS number to load from")
+        self.declareProperty(IntArrayProperty("RunNumbers", []), 'Run numbers to load')
+        self.declareProperty("Grouping", 'None', StringListValidator(['None', '2x2', '4x4']), "Group pixels")
+        self.declareProperty(WorkspaceProperty("OutputWorkspace", "",
+                                               optional=PropertyMode.Mandatory,
+                                               direction=Direction.Output),
+                             "Output Workspace")
+
+    def validateInputs(self):
+        issues = dict()
+
+        if not self.getProperty("Filename").value:
+            if (self.getProperty("IPTS").value == Property.EMPTY_INT) or len(self.getProperty("RunNumbers").value) is 0:
+                issues["Filename"] = 'Must specify either Filename or IPTS AND RunNumbers'
+
+        return issues
+
+    def PyExec(self):
+        runs = self.getProperty("Filename").value
+
+        if not runs:
+            ipts = self.getProperty("IPTS").value
+            runs = ['/HFIR/HB2C/IPTS-{}/nexus/HB2C_{}.nxs.h5'.format(ipts, run) for run in self.getProperty("RunNumbers").value]
+
+        grouping = self.getProperty("Grouping").value
+        if grouping == 'None':
+            grouping = 1
+        else:
+            grouping = 2 if grouping == '2x2' else 4
+
+        x_dim = 480*8 // grouping
+        y_dim = 512 // grouping
+
+        number_of_runs = len(runs)
+
+        data_array = np.empty((number_of_runs, x_dim, y_dim), dtype=np.float64)
+
+        s1_array = []
+        duration_array = []
+        run_number_array = []
+        monitor_count_array = []
+
+        progress = Progress(self, 0.0, 1.0, number_of_runs+3)
+
+        for n, run in enumerate(runs):
+            progress.report('Loading: '+run)
+            with h5py.File(run, 'r') as f:
+                bc = np.zeros((512*480*8),dtype=np.int64)
+                for b in range(8):
+                    bc += np.bincount(f['/entry/bank'+str(b+1)+'_events/event_id'].value,minlength=512*480*8)
+                bc = bc.reshape((480*8, 512))
+                if grouping == 2:
+                    bc = bc[::2,::2]+bc[1::2,::2]+bc[::2,1::2]+bc[1::2,1::2]
+                elif grouping == 4:
+                    bc = (bc[::4,::4]    + bc[1::4,::4]  + bc[2::4,::4]  + bc[3::4,::4]
+                          + bc[::4,1::4] + bc[1::4,1::4] + bc[2::4,1::4] + bc[3::4,1::4]
+                          + bc[::4,2::4] + bc[1::4,2::4] + bc[2::4,2::4] + bc[3::4,2::4]
+                          + bc[::4,3::4] + bc[1::4,3::4] + bc[2::4,3::4] + bc[3::4,3::4])
+                data_array[n] = bc
+                s1_array.append(f['/entry/DASlogs/HB2C:Mot:s1.RBV/average_value'].value[0])
+                duration_array.append(float(f['/entry/duration'].value[0]))
+                run_number_array.append(float(f['/entry/run_number'].value[0]))
+                monitor_count_array.append(float(f['/entry/monitor1/total_counts'].value[0]))
+
+        progress.report('Creating MDHistoWorkspace')
+        createWS_alg = self.createChildAlgorithm("CreateMDHistoWorkspace", enableLogging=False)
+        createWS_alg.setProperty("SignalInput", data_array)
+        createWS_alg.setProperty("ErrorInput", np.sqrt(data_array))
+        createWS_alg.setProperty("Dimensionality", 3)
+        createWS_alg.setProperty("Extents", '0.5,{},0.5,{},0.5,{}'.format(y_dim+0.5, x_dim+0.5, number_of_runs+0.5))
+        createWS_alg.setProperty("NumberOfBins", '{},{},{}'.format(y_dim,x_dim,number_of_runs))
+        createWS_alg.setProperty("Names", 'y,x,scanIndex')
+        createWS_alg.setProperty("Units", 'bin,bin,number')
+        createWS_alg.execute()
+        outWS = createWS_alg.getProperty("OutputWorkspace").value
+
+        progress.report('Getting IDF')
+        # Get the instrument and some logs from the first file; assume the rest are the same
+        _tmp_ws = LoadEventNexus(runs[0], MetaDataOnly=True, EnableLogging=False)
+        # The following logs should be the same for all runs
+        RemoveLogs(_tmp_ws,
+                   KeepLogs='HB2C:Mot:detz,HB2C:Mot:detz.RBV,HB2C:Mot:s2,HB2C:Mot:s2.RBV,'
+                   'HB2C:Mot:sgl,HB2C:Mot:sgl.RBV,HB2C:Mot:sgu,HB2C:Mot:sgu.RBV,'
+                   'run_title,start_time,experiment_identifier,HB2C:CS:CrystalAlign:UBMatrix',
+                   EnableLogging=False)
+
+        try:
+            ub = np.array(re.findall(r'-?\d+\.*\d*', _tmp_ws.run().getProperty('HB2C:CS:CrystalAlign:UBMatrix').value[0]),
+                          dtype=np.float).reshape(3,3)
+            sgl = np.deg2rad(_tmp_ws.run().getProperty('HB2C:Mot:sgl.RBV').value[0]) # 'HB2C:Mot:sgl.RBV,1,0,0,-1'
+            sgu = np.deg2rad(_tmp_ws.run().getProperty('HB2C:Mot:sgu.RBV').value[0]) # 'HB2C:Mot:sgu.RBV,0,0,1,-1'
+            sgl_a = np.array([[           1,            0,           0],
+                              [           0,  np.cos(sgl), np.sin(sgl)],
+                              [           0, -np.sin(sgl), np.cos(sgl)]])
+            sgu_a = np.array([[ np.cos(sgu),  np.sin(sgu),           0],
+                              [-np.sin(sgu),  np.cos(sgu),           0],
+                              [           0,            0,           1]])
+            UB = sgl_a.dot(sgu_a).dot(ub) # Apply the Goniometer tilts to the UB matrix
+            SetUB(_tmp_ws, UB=UB, EnableLogging=False)
+        except (RuntimeError, ValueError):
+            SetUB(_tmp_ws, EnableLogging=False)
+
+        if grouping > 1:
+            _tmp_group, _, _ = CreateGroupingWorkspace(InputWorkspace=_tmp_ws, EnableLogging=False)
+
+            group_number = 0
+            for x in range(0,480*8,grouping):
+                for y in range(0,512,grouping):
+                    group_number += 1
+                    for j in range(grouping):
+                        for i in range(grouping):
+                            _tmp_group.dataY(y+i+(x+j)*512)[0] = group_number
+
+            _tmp_ws = GroupDetectors(InputWorkspace=_tmp_ws, CopyGroupingFromWorkspace=_tmp_group, EnableLogging=False)
+            DeleteWorkspace(_tmp_group, EnableLogging=False)
+
+        progress.report('Adding logs')
+
+        # Hack: ConvertToMD is needed so that a deep copy of the ExperimentInfo can happen
+        # outWS.addExperimentInfo(_tmp_ws) # This doesn't work but should, when you delete `ws` `outWS` also loses it's ExperimentInfo
+        _tmp_ws = Rebin(_tmp_ws, '0,1,2', EnableLogging=False)
+        _tmp_ws = ConvertToMD(_tmp_ws, dEAnalysisMode='Elastic', EnableLogging=False, PreprocDetectorsWS='__PreprocessedDetectorsWS')
+
+        preprocWS = mtd['__PreprocessedDetectorsWS']
+        twotheta = preprocWS.column(2)
+        azimuthal = preprocWS.column(3)
+
+        outWS.copyExperimentInfos(_tmp_ws)
+        DeleteWorkspace(_tmp_ws, EnableLogging=False)
+        DeleteWorkspace('__PreprocessedDetectorsWS', EnableLogging=False)
+        # end Hack
+
+        outWS.getExperimentInfo(0).run().addProperty('s1', s1_array, True)
+        outWS.getExperimentInfo(0).run().getProperty('s1').units = 'deg'
+        outWS.getExperimentInfo(0).run().addProperty('duration', duration_array, True)
+        outWS.getExperimentInfo(0).run().getProperty('duration').units = 'second'
+        outWS.getExperimentInfo(0).run().addProperty('run_number', run_number_array, True)
+        outWS.getExperimentInfo(0).run().addProperty('monitor_count', monitor_count_array, True)
+        outWS.getExperimentInfo(0).run().addProperty('twotheta', twotheta, True)
+        outWS.getExperimentInfo(0).run().addProperty('azimuthal', azimuthal, True)
+
+        self.setProperty("OutputWorkspace", outWS)
+
+
+AlgorithmFactory.subscribe(LoadWANDSCD)
diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
index 13b8d89573b2a874abe93a41ac14e0e844daadb1..d27e83254890e80aa7e9c35f82bc7c228e144c6d 100644
--- a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
@@ -15,6 +15,7 @@ set ( TEST_PY_FILES
   CalculateSampleTransmissionTest.py
   CheckForSampleLogsTest.py
   ConjoinSpectraTest.py
+  ConvertWANDSCDtoQTest.py
   CompareSampleLogsTest.py
   ComputeCalibrationCoefVanTest.py
   ComputeIncoherentDOSTest.py
@@ -56,6 +57,7 @@ set ( TEST_PY_FILES
   LoadNMoldyn3AsciiTest.py
   LoadNMoldyn4Ascii1DTest.py
   LoadNMoldyn4AsciiTest.py
+  LoadWANDSCDTest.py
   LRPeakSelectionTest.py
   MaskAngleTest.py
   MaskBTPTest.py
diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/ConvertWANDSCDtoQTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/ConvertWANDSCDtoQTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..0789d000686149f00827e4e2e9a897c6cf460cf4
--- /dev/null
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/ConvertWANDSCDtoQTest.py
@@ -0,0 +1,138 @@
+from __future__ import absolute_import, division, print_function
+from mantid.simpleapi import ConvertWANDSCDtoQ, CreateMDHistoWorkspace, CreateSingleValuedWorkspace, SetUB, mtd
+import unittest
+import numpy as np
+
+
+class ConvertWANDSCDtoQTest(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        def gaussian(x,y,z,x0,y0,z0,ox,oy,oz,A):
+            return A*np.exp(-(x-x0)**2/(2*ox**2)-(y-y0)**2/(2*oy**2)-(z-z0)**2/(2*oz**2))
+
+        def peaks(i,j,k):
+            return gaussian(i,j,k,16,100,50,2,2,2,20)+gaussian(i,j,k,16,150,50,1,1,1,10)
+
+        S=np.fromfunction(peaks,(32,240,100))
+
+        ConvertWANDSCDtoQTest_data=CreateMDHistoWorkspace(Dimensionality=3,Extents='0.5,32.5,0.5,240.5,0.5,100.5',
+                                                          SignalInput=S.ravel('F'),ErrorInput=np.sqrt(S.ravel('F')),
+                                                          NumberOfBins='32,240,100',Names='y,x,scanIndex',Units='bin,bin,number')
+
+        ConvertWANDSCDtoQTest_dummy = CreateSingleValuedWorkspace()
+
+        ConvertWANDSCDtoQTest_data.addExperimentInfo(ConvertWANDSCDtoQTest_dummy)
+
+        ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('s1', list(np.arange(0,50,0.5)), True)
+        ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('duration', [60.]*100, True)
+        ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('monitor_count', [120000.]*100, True)
+        ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('twotheta', list(np.linspace(np.pi*2/3,0,240).repeat(32)), True)
+        ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('azimuthal', list(np.tile(np.linspace(-0.15,0.15,32),240)), True)
+
+        SetUB(ConvertWANDSCDtoQTest_data, 5,5,7,90,90,120,u=[-1,0,1],v=[1,0,1])
+
+        # Create Normalisation workspace
+        S=np.ones((32,240,1))
+        ConvertWANDSCDtoQTest_norm=CreateMDHistoWorkspace(Dimensionality=3,Extents='0.5,32.5,0.5,240.5,0.5,1.5',SignalInput=S,ErrorInput=S,
+                                                          NumberOfBins='32,240,1',Names='y,x,scanIndex',Units='bin,bin,number')
+
+        ConvertWANDSCDtoQTest_dummy2 = CreateSingleValuedWorkspace()
+        ConvertWANDSCDtoQTest_norm.addExperimentInfo(ConvertWANDSCDtoQTest_dummy2)
+        ConvertWANDSCDtoQTest_norm.getExperimentInfo(0).run().addProperty('monitor_count', [100000.], True)
+
+    @classmethod
+    def tearDownClass(cls):
+        [mtd.remove(ws) for ws in ['ConvertWANDSCDtoQTest_data',
+                                   'ConvertWANDSCDtoQTest_dummy'
+                                   'ConvertWANDSCDtoQTest_norm',
+                                   'ConvertWANDSCDtoQTest_dummy2']]
+
+    def test_Q(self):
+        ConvertWANDSCDtoQTest_out=ConvertWANDSCDtoQ('ConvertWANDSCDtoQTest_data',BinningDim0='-8.08,8.08,101',
+                                                    BinningDim1='-0.88,0.88,11',BinningDim2='-8.08,8.08,101',NormaliseBy='None')
+
+        self.assertTrue(ConvertWANDSCDtoQTest_out)
+
+        s = ConvertWANDSCDtoQTest_out.getSignalArray()
+        self.assertAlmostEqual(np.nanmax(s), 8.97233331213612)
+        self.assertAlmostEqual(np.nanargmax(s), 22780)
+
+        self.assertEquals(ConvertWANDSCDtoQTest_out.getNumDims(), 3)
+        self.assertEquals(ConvertWANDSCDtoQTest_out.getNPoints(), 112211)
+
+        d0 = ConvertWANDSCDtoQTest_out.getDimension(0)
+        self.assertEquals(d0.name, 'Q_sample_x')
+        self.assertEquals(d0.getNBins(), 101)
+        self.assertAlmostEquals(d0.getMinimum(), -8.08, 5)
+        self.assertAlmostEquals(d0.getMaximum(), 8.08, 5)
+
+        d1 = ConvertWANDSCDtoQTest_out.getDimension(1)
+        self.assertEquals(d1.name, 'Q_sample_y')
+        self.assertEquals(d1.getNBins(), 11)
+        self.assertAlmostEquals(d1.getMinimum(), -0.88, 5)
+        self.assertAlmostEquals(d1.getMaximum(), 0.88, 5)
+
+        d2 = ConvertWANDSCDtoQTest_out.getDimension(2)
+        self.assertEquals(d2.name, 'Q_sample_z')
+        self.assertEquals(d2.getNBins(), 101)
+        self.assertAlmostEquals(d2.getMinimum(), -8.08, 5)
+        self.assertAlmostEquals(d2.getMaximum(), 8.08, 5)
+
+        self.assertEqual(ConvertWANDSCDtoQTest_out.getNumExperimentInfo(), 1)
+
+        ConvertWANDSCDtoQTest_out.delete()
+
+    def test_Q_norm(self):
+        ConvertWANDSCDtoQTest_out = ConvertWANDSCDtoQ('ConvertWANDSCDtoQTest_data',NormalisationWorkspace='ConvertWANDSCDtoQTest_norm',
+                                                      BinningDim0='-8.08,8.08,101',BinningDim1='-0.88,0.88,11',BinningDim2='-8.08,8.08,101')
+
+        s = ConvertWANDSCDtoQTest_out.getSignalArray()
+        self.assertAlmostEqual(np.nanmax(s), 7.476944426780101)
+        self.assertAlmostEqual(np.nanargmax(s), 22780)
+
+        ConvertWANDSCDtoQTest_out.delete()
+
+    def test_HKL_norm_and_KeepTemporary(self):
+        ConvertWANDSCDtoQTest_out = ConvertWANDSCDtoQ('ConvertWANDSCDtoQTest_data',NormalisationWorkspace='ConvertWANDSCDtoQTest_norm',
+                                                      Frame='HKL',KeepTemporaryWorkspaces=True,BinningDim0='-8.08,8.08,101',
+                                                      BinningDim1='-8.08,8.08,101',BinningDim2='-8.08,8.08,101',
+                                                      Uproj='1,1,0',Vproj='1,-1,0',Wproj='0,0,1')
+
+        self.assertTrue(ConvertWANDSCDtoQTest_out)
+        self.assertTrue(mtd.doesExist('ConvertWANDSCDtoQTest_out'))
+        self.assertTrue(mtd.doesExist('ConvertWANDSCDtoQTest_out_data'))
+        self.assertTrue(mtd.doesExist('ConvertWANDSCDtoQTest_out_normalization'))
+
+        s = ConvertWANDSCDtoQTest_out.getSignalArray()
+        self.assertAlmostEqual(np.nanmax(s), 4.646855396509936)
+        self.assertAlmostEqual(np.nanargmax(s), 443011)
+
+        self.assertEquals(ConvertWANDSCDtoQTest_out.getNumDims(), 3)
+        self.assertEquals(ConvertWANDSCDtoQTest_out.getNPoints(), 101**3)
+
+        d0 = ConvertWANDSCDtoQTest_out.getDimension(0)
+        self.assertEquals(d0.name, '[H,H,0]')
+        self.assertEquals(d0.getNBins(), 101)
+        self.assertAlmostEquals(d0.getMinimum(), -8.08, 5)
+        self.assertAlmostEquals(d0.getMaximum(), 8.08, 5)
+
+        d1 = ConvertWANDSCDtoQTest_out.getDimension(1)
+        self.assertEquals(d1.name, '[H,-H,0]')
+        self.assertEquals(d1.getNBins(), 101)
+        self.assertAlmostEquals(d1.getMinimum(), -8.08, 5)
+        self.assertAlmostEquals(d1.getMaximum(), 8.08, 5)
+
+        d2 = ConvertWANDSCDtoQTest_out.getDimension(2)
+        self.assertEquals(d2.name, '[0,0,L]')
+        self.assertEquals(d2.getNBins(), 101)
+        self.assertAlmostEquals(d2.getMinimum(), -8.08, 5)
+        self.assertAlmostEquals(d2.getMaximum(), 8.08, 5)
+
+        self.assertEqual(ConvertWANDSCDtoQTest_out.getNumExperimentInfo(), 1)
+
+        ConvertWANDSCDtoQTest_out.delete()
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/LoadWANDSCDTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/LoadWANDSCDTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..616941e7a415f69f5db82e89722d22ce6a990972
--- /dev/null
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/LoadWANDSCDTest.py
@@ -0,0 +1,58 @@
+from __future__ import absolute_import, division, print_function
+from mantid.simpleapi import LoadWANDSCD
+import unittest
+
+
+class LoadWANDTest(unittest.TestCase):
+
+    def test(self):
+        LoadWANDTest_ws = LoadWANDSCD('HB2C_7000.nxs.h5,HB2C_7001.nxs.h5')
+        self.assertTrue(LoadWANDTest_ws)
+        self.assertEquals(LoadWANDTest_ws.getNumDims(), 3)
+        self.assertEquals(LoadWANDTest_ws.getNPoints(), 1966080*2)
+        self.assertEquals(LoadWANDTest_ws.getSignalArray().max(), 7)
+
+        d0 = LoadWANDTest_ws.getDimension(0)
+        self.assertEquals(d0.name, 'y')
+        self.assertEquals(d0.getNBins(), 512)
+        self.assertEquals(d0.getMinimum(), 0.5)
+        self.assertEquals(d0.getMaximum(), 512.5)
+
+        d1 = LoadWANDTest_ws.getDimension(1)
+        self.assertEquals(d1.name, 'x')
+        self.assertEquals(d1.getNBins(), 3840)
+        self.assertEquals(d1.getMinimum(), 0.5)
+        self.assertEquals(d1.getMaximum(), 3840.5)
+
+        d2 = LoadWANDTest_ws.getDimension(2)
+        self.assertEquals(d2.name, 'scanIndex')
+        self.assertEquals(d2.getNBins(), 2)
+        self.assertEquals(d2.getMinimum(), 0.5)
+        self.assertEquals(d2.getMaximum(), 2.5)
+
+        self.assertEqual(LoadWANDTest_ws.getNumExperimentInfo(), 1)
+        self.assertEquals(LoadWANDTest_ws.getExperimentInfo(0).getInstrument().getName(), 'WAND')
+
+        run = LoadWANDTest_ws.getExperimentInfo(0).run()
+        s1 = run.getProperty('s1').value
+        self.assertEqual(len(s1), 2)
+        self.assertEqual(s1[0], -142.6)
+        self.assertEqual(s1[1], -142.5)
+        run_number = run.getProperty('run_number').value
+        self.assertEqual(len(run_number), 2)
+        self.assertEqual(run_number[0], 7000)
+        self.assertEqual(run_number[1], 7001)
+        monitor_count=run.getProperty('monitor_count').value
+        self.assertEqual(len(monitor_count), 2)
+        self.assertEqual(monitor_count[0], 907880)
+        self.assertEqual(monitor_count[1], 908651)
+        duration = run.getProperty('duration').value
+        self.assertEqual(len(duration), 2)
+        self.assertAlmostEqual(duration[0], 40.05, 5)
+        self.assertAlmostEqual(duration[1], 40.05, 5)
+
+        LoadWANDTest_ws.delete()
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/Testing/Data/SystemTest/HB2C_WANDSCD_data.nxs.md5 b/Testing/Data/SystemTest/HB2C_WANDSCD_data.nxs.md5
new file mode 100644
index 0000000000000000000000000000000000000000..ca989f82037f45c75bdb41b7b162b9218b271380
--- /dev/null
+++ b/Testing/Data/SystemTest/HB2C_WANDSCD_data.nxs.md5
@@ -0,0 +1 @@
+486a2e500a75b47d857cb9c8ef7cb973
diff --git a/Testing/Data/SystemTest/HB2C_WANDSCD_norm.nxs.md5 b/Testing/Data/SystemTest/HB2C_WANDSCD_norm.nxs.md5
new file mode 100644
index 0000000000000000000000000000000000000000..b3a923cf4d363e28266c81775af3643a16c04363
--- /dev/null
+++ b/Testing/Data/SystemTest/HB2C_WANDSCD_norm.nxs.md5
@@ -0,0 +1 @@
+4e902cca48f73de3414490073d2c73b5
diff --git a/Testing/Data/UnitTest/HB2C_7001.nxs.h5.md5 b/Testing/Data/UnitTest/HB2C_7001.nxs.h5.md5
new file mode 100644
index 0000000000000000000000000000000000000000..5249348146fd437a6a53d3c09018b5a5306767c0
--- /dev/null
+++ b/Testing/Data/UnitTest/HB2C_7001.nxs.h5.md5
@@ -0,0 +1 @@
+1faee25c35ae34bcdc13618224547a23
diff --git a/Testing/SystemTests/tests/analysis/ConvertWANDSCDtoQTest.py b/Testing/SystemTests/tests/analysis/ConvertWANDSCDtoQTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..fc1adb01fa57034410725c6d19d676fef3714702
--- /dev/null
+++ b/Testing/SystemTests/tests/analysis/ConvertWANDSCDtoQTest.py
@@ -0,0 +1,52 @@
+import stresstesting
+import numpy as np
+from mantid.simpleapi import *
+
+
+class ConvertWANDSCDtoQTest(stresstesting.MantidStressTest):
+    def requiredMemoryMB(self):
+            return 8000
+
+    def runTest(self):
+        LoadMD('HB2C_WANDSCD_data.nxs', OutputWorkspace='ConvertWANDSCDtoQTest_data')
+        LoadMD('HB2C_WANDSCD_norm.nxs', OutputWorkspace='ConvertWANDSCDtoQTest_norm')
+        ConvertWANDSCDtoQTest_Q=ConvertWANDSCDtoQ(InputWorkspace='ConvertWANDSCDtoQTest_data',
+                                                  NormalisationWorkspace='ConvertWANDSCDtoQTest_norm')
+
+        ConvertWANDSCDtoQTest_peaks=FindPeaksMD(InputWorkspace=ConvertWANDSCDtoQTest_Q, PeakDistanceThreshold=2,
+                                                CalculateGoniometerForCW=True, Wavelength=1.488)
+        FindUBUsingLatticeParameters(ConvertWANDSCDtoQTest_peaks, a=5.64, b=5.64, c=5.64, alpha=90, beta=90, gamma=90)
+        UB = ConvertWANDSCDtoQTest_peaks.sample().getOrientedLattice().getUB()
+        self.assertTrue(np.allclose(UB, [[-2.73152432e-17,  1.77061974e-01, -9.27942487e-03],
+                                         [ 1.77304965e-01,  0.00000000e+00,  0.00000000e+00],
+                                         [ 1.23032284e-17, -9.27942487e-03, -1.77061974e-01]]))
+
+        ConvertWANDSCDtoQ(InputWorkspace='ConvertWANDSCDtoQTest_data',
+                          NormalisationWorkspace='ConvertWANDSCDtoQTest_norm',
+                          UBWorkspace=ConvertWANDSCDtoQTest_peaks,
+                          Frame='HKL',
+                          BinningDim0='-0.62,0.62,31',
+                          BinningDim1='-2.02,7.02,226',
+                          BinningDim2='-6.52,2.52,226',
+                          OutputWorkspace='ConvertWANDSCDtoQTest_HKL')
+
+    def validate(self):
+        results = 'ConvertWANDSCDtoQTest_HKL'
+        reference = 'ConvertWANDSCDtoQTest_HKL.nxs'
+
+        Load(Filename=reference,OutputWorkspace=reference)
+
+        checker = AlgorithmManager.create("CompareMDWorkspaces")
+        checker.setLogging(True)
+        checker.setPropertyValue("Workspace1",results)
+        checker.setPropertyValue("Workspace2",reference)
+        checker.setPropertyValue("Tolerance", "1e-7")
+
+        checker.execute()
+        if checker.getPropertyValue("Equals") != "1":
+            print(" Workspaces do not match, result: ",checker.getPropertyValue("Result"))
+            print(self.__class__.__name__)
+            SaveMD(InputWorkspace=results,Filename=self.__class__.__name__+'-mismatch.nxs')
+            return False
+
+        return True
diff --git a/Testing/SystemTests/tests/analysis/reference/ConvertWANDSCDtoQTest_HKL.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/ConvertWANDSCDtoQTest_HKL.nxs.md5
new file mode 100644
index 0000000000000000000000000000000000000000..722c5ac489e20912b8e27d3b6fb16ecc7a0a34e5
--- /dev/null
+++ b/Testing/SystemTests/tests/analysis/reference/ConvertWANDSCDtoQTest_HKL.nxs.md5
@@ -0,0 +1 @@
+8a827de4e0d686c66cc2f1593d4745fd
diff --git a/docs/source/algorithms/ConvertWANDSCDtoQ-v1.rst b/docs/source/algorithms/ConvertWANDSCDtoQ-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..8dcd25ac2f9f1647699ced7b764a688d088448b2
--- /dev/null
+++ b/docs/source/algorithms/ConvertWANDSCDtoQ-v1.rst
@@ -0,0 +1,94 @@
+.. algorithm::
+
+.. summary::
+
+.. relatedalgorithms::
+
+.. properties::
+
+Description
+-----------
+
+This algorithm will convert the output of :ref:`algm-LoadWANDSCD` in
+either Q or HKL space. :ref:`algm-FindPeaksMD` can be run on the
+output Q sample space, then the UB can be found and used to then
+convert to HKL. The default binning ranges are good for converting to
+Q sample with the default wavelength.
+
+The normalization is calculated in the same way as
+:ref:`algm-MDNormSCD` but with the solid angle and flux coming from
+the NormalisationWorkspace, normally vanadium. A brief introduction to
+the multi-dimensional data normalization can be found :ref:`here
+<MDNorm>`.
+
+When converting to HKL it will use the UB matrix from the UBWorkspace
+if provided otherwise it will use the UB matrix from the
+InputWorkspace. Uproj, Vproj and Wproj are only used when converting
+to HKL Frame.
+
+If the KeepTemporaryWorkspaces option is True the data and the
+normalization in addition to the nomalized data will be
+outputted. This allows you to run separate instances of
+ConvertWANDSCDtoQ and combine the results. They will have names
+"ws_data" and "ws_normalization" respectively.
+
+Usage
+-----
+
+**Convert to Q**
+
+.. code-block:: python
+
+    # Load Data and normalisation
+    LoadWANDSCD(IPTS=7776, RunNumbers=26509, OutputWorkspace='norm',Grouping='4x4') # Vanadium
+    LoadWANDSCD(IPTS=7776, RunNumbers='26640-27944', OutputWorkspace='data',Grouping='4x4')
+    ConvertWANDSCDtoQ(InputWorkspace='data',
+                      NormalisationWorkspace='norm',
+                      OutputWorkspace='Q',
+                      BinningDim1='-1,1,1')
+
+    # Plot workspace
+    import matplotlib.pyplot as plt
+    from mantid import plots
+    fig, ax = plt.subplots(subplot_kw={'projection':'mantid'})
+    c = ax.pcolormesh(mtd['Q'], vmax=1)
+    cbar=fig.colorbar(c)
+    cbar.set_label('Intensity (arb. units)')
+    #fig.savefig('ConvertWANDSCDtoQ_Q.png')
+
+Output:
+
+.. figure:: /images/ConvertWANDSCDtoQ_Q.png
+
+**Convert to HKL**
+
+.. code-block:: python
+
+    # Load Data and normalisation
+    LoadWANDSCD(IPTS=7776, RunNumbers=26509, OutputWorkspace='norm',Grouping='4x4') # Vanadium
+    LoadWANDSCD(IPTS=7776, RunNumbers='26640-27944', OutputWorkspace='data',Grouping='4x4')
+    SetUB('data', UB='0,0.1770619741,-0.00927942487,0.177304965,0,0,0,-0.00927942487,-0.177061974')
+    ConvertWANDSCDtoQ(InputWorkspace='data',
+                      NormalisationWorkspace='norm',
+                      OutputWorkspace='HKL',
+                      Frame='HKL',
+                      BinningDim0='-1,1,1',
+                      BinningDim1='-2.02,7.02,226',
+                      BinningDim2='-6.52,2.52,226')
+
+    # Plot workspace
+    import matplotlib.pyplot as plt
+    from mantid import plots
+    fig, ax = plt.subplots(subplot_kw={'projection':'mantid'})
+    c = ax.pcolormesh(mtd['HKL'], vmax=1)
+    cbar=fig.colorbar(c)
+    cbar.set_label('Intensity (arb. units)')
+    #fig.savefig('ConvertWANDSCDtoQ_HKL.png')
+
+Output:
+
+.. figure:: /images/ConvertWANDSCDtoQ_HKL.png
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/algorithms/LoadWANDSCD-v1.rst b/docs/source/algorithms/LoadWANDSCD-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..5da397a27378291c5217ab607cff71372d9ea311
--- /dev/null
+++ b/docs/source/algorithms/LoadWANDSCD-v1.rst
@@ -0,0 +1,143 @@
+.. algorithm::
+
+.. summary::
+
+.. relatedalgorithms::
+
+.. properties::
+
+Description
+-----------
+
+This algorithm will load a series of run into a MDHistoWorkspace that
+has dimensions x and y detector pixels vs scanIndex. The scanIndex is
+the omega rotation of the sample. The instrument from the first run
+only will be copied to the OutputWorkspace. In addition the s1 (omega
+rotation), duration, run_number and monitor count is read from every
+file and included in the logs of the OutputWorkspace.
+
+If the "HB2C:CS:CrystalAlign:UBMatrix" property exist it will be
+converted into the OrientedLattice on the OutputWorkspace. The
+goniometer tilts (sgu and sgl) are combined into the UB Matrix so that
+only omega (s1) needs to be taken into account during rotation.
+
+This algorithm doesn't use Mantid loaders but instead h5py and numpy
+to load and integrate the events.
+
+There is a grouping option to group pixels by either 2x2 or 4x4 which
+will help in reducing memory usage and speed up the later reduction
+steps. In most cases you will not see a difference in reduced data
+with 4x4 pixel grouping.
+
+The loaded workspace is designed to be the input to
+:ref:`algm-ConvertWANDSCDtoQ`.
+
+Usage
+-----
+
+**Load one file, Vanadium for normalisation**
+
+.. code-block:: python
+
+    norm = LoadWANDSCD(IPTS=7776, RunNumbers=26509)
+    print(repr(norm))
+
+Output:
+
+.. code-block:: none
+
+    MDHistoWorkspace
+    Title:
+    Dim 0: (y) 0.5 to 512.5 in 512 bins
+    Dim 1: (x) 0.5 to 3840.5 in 3840 bins
+    Dim 2: (scanIndex) 0.5 to 1.5 in 1 bins
+
+    Inelastic: ki-kf
+    Instrument: ...
+
+    Run start: 2018-Mar-12 17:10:59
+    Run end:  not available
+    Sample: a 1.0, b 1.0, c 1.0; alpha 90, beta 90, gamma 90
+
+
+**Load multiple data file**
+
+.. code-block:: python
+
+    data = LoadWANDSCD(IPTS=7776, RunNumbers='26640-27944')
+    print("Memory used: {}GiB".format(data.getMemorySize()/2**30))
+    print(repr(data))
+    print('s1 = {}'.format(data.getExperimentInfo(0).run().getProperty('s1').value[0:10]))
+    print('monitor_counts = {}'.format(data.getExperimentInfo(0).run().getProperty('monitor_counts').value[0:10]))
+    print('duration = {}'.format(data.getExperimentInfo(0).run().getProperty('duration').value[0:10]))
+    print('run_number = {}'.format(data.getExperimentInfo(0).run().getProperty('run_number').value[0:10]))
+
+Output:
+
+.. code-block:: none
+
+    Memory used: 59GB
+
+    MDHistoWorkspace
+    Title:
+    Dim 0: (y) 0.5 to 512.5 in 512 bins
+    Dim 1: (x) 0.5 to 3840.5 in 3840 bins
+    Dim 2: (scanIndex) 0.5 to 1305.5 in 1305 bins
+    Inelastic: ki-kf
+    Instrument: ...
+    Run start: 2018-May-02 13:34:10
+    Run end:  not available
+    Sample: a 5.7, b 5.7, c 5.6; alpha 93, beta 90, gamma 98
+
+    s2 = [-180,-179.9,-179.8,-179.7,-179.6,-179.5,-179.4,-179.3,-179.2,-179.1]
+    monitor_count = [44571,44598,44567,44869,44453,44238,44611,44120,44762,44658]
+    duration = [2.05,2.05,2.03333,2.05,2.03333,2.03333,2.05,2.01667,2.05,2.05]
+    run_number = [26640,26641,26642,26643,26644,26645,26646,26647,26648,26649]
+
+
+**Load with different grouping comparing memory usage**
+
+.. code-block:: python
+
+    data = LoadWANDSCD(IPTS=7776, RunNumbers='26640-27944')
+    data_2x2 = LoadWANDSCD(IPTS=7776, RunNumbers='26640-27944', Grouping='2x2')
+    data_4x4 = LoadWANDSCD(IPTS=7776, RunNumbers='26640-27944', Grouping='4x4')
+    print("Memory used by {}: {}GiB".format(data,data.getMemorySize()/2**30))
+    print("Memory used by {}: {}GiB".format(data_2x2,data_2x2.getMemorySize()/2**30))
+    print("Memory used by {}: {}GiB".format(data_4x4,data_4x4.getMemorySize()/2**30))
+    print(repr(data_4x4))
+
+    # Integrate y and plot
+    data_integrated = IntegrateMDHistoWorkspace('data_4x4', P1Bin='0,129')
+    import matplotlib.pyplot as plt
+    from mantid import plots
+    fig, ax = plt.subplots(subplot_kw={'projection':'mantid'})
+    c = ax.pcolormesh(data_integrated, vmax=100)
+    cbar=fig.colorbar(c)
+    cbar.set_label('Intensity (arb. units)')
+    #fig.savefig('LoadWANDSCD.png')
+
+Output:
+
+.. code-block:: none
+
+    Memory used by data: 59GiB
+    Memory used by data_2x2: 14GiB
+    Memory used by data_4x4: 3GiB
+
+    MDHistoWorkspace
+    Title:
+    Dim 0: (y) 0.5 to 128.5 in 128 bins
+    Dim 1: (x) 0.5 to 960.5 in 960 bins
+    Dim 2: (scanIndex) 0.5 to 1305.5 in 1305 bins
+    Inelastic: ki-kf
+    Instrument: ...
+    Run start: 2018-May-02 13:34:10
+    Run end:  not available
+    Sample: a 5.7, b 5.7, c 5.6; alpha 93, beta 90, gamma 98
+
+.. figure:: /images/LoadWANDSCD.png
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/images/ConvertWANDSCDtoQ_HKL.png b/docs/source/images/ConvertWANDSCDtoQ_HKL.png
new file mode 100644
index 0000000000000000000000000000000000000000..f65b355198c3cbeca851123bc6990fa7b41a70c0
Binary files /dev/null and b/docs/source/images/ConvertWANDSCDtoQ_HKL.png differ
diff --git a/docs/source/images/ConvertWANDSCDtoQ_Q.png b/docs/source/images/ConvertWANDSCDtoQ_Q.png
new file mode 100644
index 0000000000000000000000000000000000000000..bcc41290cc3d5fb08f7b9a7007ba8a7c1b4030e2
Binary files /dev/null and b/docs/source/images/ConvertWANDSCDtoQ_Q.png differ
diff --git a/docs/source/images/LoadWANDSCD.png b/docs/source/images/LoadWANDSCD.png
new file mode 100644
index 0000000000000000000000000000000000000000..168452904a3b3ee8988c1452fda7bd32702bfbe0
Binary files /dev/null and b/docs/source/images/LoadWANDSCD.png differ
diff --git a/docs/source/release/v3.13.0/diffraction.rst b/docs/source/release/v3.13.0/diffraction.rst
index 5195e091a2593f3d1b298c7b5dd84d438de5c001..024150a0f67358f9e6739dcab963e3fe2e782fef 100644
--- a/docs/source/release/v3.13.0/diffraction.rst
+++ b/docs/source/release/v3.13.0/diffraction.rst
@@ -61,6 +61,8 @@ Single Crystal Diffraction
 
 - :ref:`CreatePeaksWorkspace <algm-CreatePeaksWorkspace>` now accepts MD workspaces as input.
 
+- New algorithms :ref:`LoadWANDSCD <algm-LoadWANDSCD>` and :ref:`ConvertWANDSCDtoQ <algm-ConvertWANDSCDtoQ>` to load single crystal HB2C data and convert it to either Q-sample or HKL space.
+
 Improvements
 ############