diff --git a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/Quest.py b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/Quest.py
deleted file mode 100644
index 5f15cb3c2fb160048a520cf1a8bbd38c07254126..0000000000000000000000000000000000000000
--- a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/Quest.py
+++ /dev/null
@@ -1,107 +0,0 @@
-#pylint: disable=no-init
-from mantid.api import PythonAlgorithm, AlgorithmFactory
-from mantid.kernel import StringListValidator, StringMandatoryValidator
-from mantid.simpleapi import *
-from mantid import config, logger, mtd
-import os
-
-class Quest(PythonAlgorithm):
-
-    def category(self):
-        return "Workflow\\MIDAS"
-
-    def summary(self):
-        return "This is a variation of the stretched exponential option of Quasi."
-
-    def PyInit(self):
-        self.declareProperty(name='InputType',defaultValue='File',validator=StringListValidator(['File','Workspace']),
-                             doc='Origin of data input - File (.nxs) or Workspace')
-        self.declareProperty(name='Instrument',defaultValue='iris',validator=StringListValidator(['irs','iris','osi','osiris']),
-                             doc='Instrument')
-        self.declareProperty(name='Analyser',defaultValue='graphite002',validator=StringListValidator(['graphite002','graphite004']),
-                             doc='Analyser & reflection')
-        self.declareProperty(name='SamNumber',defaultValue='',validator=StringMandatoryValidator(), doc='Sample run number')
-        self.declareProperty(name='ResInputType',defaultValue='File',validator=StringListValidator(['File','Workspace']),
-                             doc='Origin of res input - File (_res.nxs) or Workspace')
-        self.declareProperty(name='ResNumber',defaultValue='',validator=StringMandatoryValidator(), doc='Resolution run number')
-        self.declareProperty(name='ResNormInputType',defaultValue='File',validator=StringListValidator(['File','Workspace']),
-                             doc='Origin of ResNorm input - File (_red.nxs) or Workspace')
-        self.declareProperty(name='ResNormNumber',defaultValue='',validator=StringMandatoryValidator(), doc='ResNorm run number')
-        self.declareProperty(name='ElasticOption',defaultValue=True, doc='Include elastic peak in fit')
-        self.declareProperty(name='BackgroundOption',defaultValue='Sloping',validator=StringListValidator(['Sloping','Flat','Zero']),
-                             doc='Form of background to fit')
-        self.declareProperty(name='EnergyMin', defaultValue=-0.5, doc='Minimum energy for fit. Default=-0.5')
-        self.declareProperty(name='EnergyMax', defaultValue=0.5, doc='Maximum energy for fit. Default=0.5')
-        self.declareProperty(name='SamBinning', defaultValue=1, doc='Binning value(integer) for sample. Default=1')
-        self.declareProperty(name='NumberSigma', defaultValue=50, doc='Number of sigma values. Default=50')
-        self.declareProperty(name='NumberBeta', defaultValue=30, doc='Number of beta values. Default=30')
-        self.declareProperty(name='Sequence',defaultValue=True, doc='Switch Sequence Off/On')
-        self.declareProperty(name='Plot',defaultValue='None',validator=StringListValidator(['None','Sigma','Beta','All']),
-                             doc='Plot options')
-        self.declareProperty(name='Verbose',defaultValue=True, doc='Switch Verbose Off/On')
-        self.declareProperty(name='Save',defaultValue=False, doc='Switch Save result to nxs file Off/On')
-
-    #pylint: disable=too-many-locals,too-many-branches
-    def PyExec(self):
-        from IndirectImport import run_f2py_compatibility_test, is_supported_f2py_platform
-
-        if is_supported_f2py_platform():
-            import IndirectBayes as Main
-
-        run_f2py_compatibility_test()
-
-        self.log().information('Quest input')
-        inType = self.getPropertyValue('InputType')
-        prefix = self.getPropertyValue('Instrument')
-        ana = self.getPropertyValue('Analyser')
-        sam = self.getPropertyValue('SamNumber')
-        rinType = self.getPropertyValue('ResInputType')
-        res = self.getPropertyValue('ResNumber')
-        elastic = self.getProperty('ElasticOption').value
-        bgd = self.getPropertyValue('BackgroundOption')
-        emin = self.getPropertyValue('EnergyMin')
-        emax = self.getPropertyValue('EnergyMax')
-        nbin = self.getPropertyValue('SamBinning')
-        nbins = [nbin, 1]
-        nbet = self.getProperty('NumberBeta').value
-        nsig = self.getProperty('NumberSigma').value
-        nbs = [nbet, nsig]
-
-        sname = prefix+sam+'_'+ana + '_red'
-        rname = prefix+res+'_'+ana + '_res'
-        erange = [float(emin), float(emax)]
-        if elastic:
-            o_el = 1
-        else:
-            o_el = 0
-        fitOp = [o_el, bgd, 0, 0]
-        loopOp = self.getProperty('Sequence').value
-        verbOp = self.getProperty('Verbose').value
-        plotOp = self.getPropertyValue('Plot')
-        saveOp = self.getProperty('Save').value
-
-        workdir = config['defaultsave.directory']
-        if not os.path.isdir(workdir):
-            workdir = os.getcwd()
-            logger.information('Default Save directory is not set. Defaulting to current working Directory: ' + workdir)
-
-        if inType == 'File':
-            spath = os.path.join(workdir, sname+'.nxs')        # path name for sample nxs file
-            LoadNexusProcessed(Filename=spath, OutputWorkspace=sname)
-            Smessage = 'Sample from File : '+spath
-        else:
-            Smessage = 'Sample from Workspace : '+sname
-
-        if rinType == 'File':
-            rpath = os.path.join(workdir, rname+'.nxs')        # path name for res nxs file
-            LoadNexusProcessed(Filename=rpath, OutputWorkspace=rname)
-            Rmessage = 'Resolution from File : '+rpath
-        else:
-            Rmessage = 'Resolution from Workspace : '+rname
-
-        if verbOp:
-            logger.notice(Smessage)
-            logger.notice(Rmessage)
-        Main.QuestRun(sname,rname,nbs,erange,nbins,fitOp,loopOp,plotOp,saveOp)
-
-AlgorithmFactory.subscribe(Quest)         # Register algorithm with Mantid
diff --git a/Testing/SystemTests/tests/analysis/ISISIndirectBayesTest.py b/Testing/SystemTests/tests/analysis/ISISIndirectBayesTest.py
index b2cf8faa7b541b6c4a3185b4f61a6be6349af43d..5cede9caa98d97d840d67bc1e04e30c1a8bed851 100644
--- a/Testing/SystemTests/tests/analysis/ISISIndirectBayesTest.py
+++ b/Testing/SystemTests/tests/analysis/ISISIndirectBayesTest.py
@@ -101,27 +101,30 @@ class QuestTest(stresstesting.MantidStressTest):
         return not platform.system() == "Windows"
 
     def runTest(self):
-        import IndirectBayes as Main
-        nbins = [1, 1]
-        nbs = [50, 30]
         sname = 'irs26176_graphite002_red'
         rname = 'irs26173_graphite002_res'
-        erange = [-0.5, 0.5]
-        fitOp = [True, 'Sloping', False, False] #elastic, background, width, resnorm
-        loopOp = False
-        plotOp = 'None'
-        saveOp = False
 
         spath = sname+'.nxs'   # path name for sample nxs file
-        LoadNexusProcessed(Filename=spath, OutputWorkspace=sname)
+        sample = LoadNexusProcessed(Filename=spath, OutputWorkspace=sname)
         rpath = rname+'.nxs'    # path name for res nxs file
-        LoadNexusProcessed(Filename=rpath, OutputWorkspace=rname)
-        Main.QuestRun(sname,rname,nbs,erange,nbins,fitOp,loopOp,plotOp,saveOp)
+        res = LoadNexusProcessed(Filename=rpath, OutputWorkspace=rname)
+        BayesStretch(SampleWorkspace=sample,
+                     ResolutionWorkspace=res,
+                     EMin=-0.5,
+                     EMax=0.5,
+                     SampleBins=1,
+                     Elastic=True,
+                     Background='Sloping',
+                     NumberSigma=30,
+                     NumberBeta=50,
+                     Loop=False,
+                     OutputWorkspaceFit='fit_group',
+                     OutputWorkspaceContour='contour_group')
 
     def validate(self):
         self.tolerance = 1e-1
         self.disableChecking.append('SpectraMap')
-        return 'irs26176_graphite002_Qst_Fit','ISISIndirectBayes_QuestTest.nxs'
+        return 'irs26176_graphite002_Stretch_Fit','ISISIndirectBayes_QuestTest.nxs'
 
     def cleanup(self):
         filenames = ['irs26176_graphite002_Qst.lpt','irs26176_graphite002_Qss.ql2',
diff --git a/docs/source/algorithms/Quest-v1.rst b/docs/source/algorithms/Quest-v1.rst
deleted file mode 100644
index 7ffd94f7162524b0f3f5ad6b4f9c32cce158ae7e..0000000000000000000000000000000000000000
--- a/docs/source/algorithms/Quest-v1.rst
+++ /dev/null
@@ -1,35 +0,0 @@
-.. algorithm::
-
-.. summary::
-
-.. alias::
-
-.. properties::
-
-Description
------------
-
-This is a variation of the stretched exponential option of
-`Quasi <http://www.mantidproject.org/IndirectBayes:Quasi>`__. For each spectrum a fit is performed
-for a grid of :math:`\beta` and :math:`\sigma` values. The distribution of goodness of fit values
-is plotted.
-
-This routine was originally part of the MODES package. Note that this algorithm
-uses F2Py and is currently only supported on windows.
-
-Usage
------
-**Example - a basic example using Quest to fit a reduced workspace.**
-
-.. code-block:: python
-
-   sam = LoadNexusProcessed("irs26173_graphite002_red.nxs", OutputWorkspace='irs26173_graphite002_red')
-   res = LoadNexusProcessed("irs26173_graphite002_res.nxs", OutputWorkspace='irs26173_graphite002_res')
-   van = LoadNexusProcessed("irs26176_graphite002_red.nxs", OutputWorkspace='irs26176_graphite002_red')
-
-   ResNorm(VanNumber='26176', ResNumber='26173', InputType='File', ResInputType='File', Instrument='irs', Analyser='graphite002', Plot='None', Version=1)
-   Quest(SamNumber='26176', ResNumber='26173', ResNormNumber='26176', InputType='File', ResInputType='File', ResNormInputType='Workspace', Instrument='irs', Analyser='graphite002')
-
-.. categories::
-
-.. sourcelink::
diff --git a/scripts/Inelastic/IndirectBayes.py b/scripts/Inelastic/IndirectBayes.py
index 4f0bf3c7135f580476e502c960c28aca1419a821..105bb67b4e9583b7e754a1ee5c906663acf5db29 100644
--- a/scripts/Inelastic/IndirectBayes.py
+++ b/scripts/Inelastic/IndirectBayes.py
@@ -12,7 +12,6 @@ if is_supported_f2py_platform():
     QLr     = import_f2py("QLres")
     QLd     = import_f2py("QLdata")
     Qse     = import_f2py("QLse")
-    Que     = import_f2py("Quest")
     resnorm = import_f2py("ResNorm")
 else:
     unsupported_message()
@@ -68,203 +67,6 @@ def GetXYE(inWS,n,array_len):
     E=PadArray(Ein,array_len)
     return N,X,Y,E
 
-# Quest programs
-def CheckBetSig(nbs):
-    Nsig = int(nbs[1])
-    if Nsig == 0:
-        raise ValueError('Number of sigma points is Zero')
-    if Nsig > 200:
-        raise ValueError('Max number of sigma points is 200')
-
-    Nbet = int(nbs[0])
-    if Nbet == 0:
-        raise ValueError('Number of beta points is Zero')
-    if Nbet > 200:
-        raise ValueError('Max number of beta points is 200')
-
-    return Nbet,Nsig
-
-def QuestRun(samWS,resWS,nbs,erange,nbins,Fit,Loop,Plot,Save):
-    StartTime('Quest')
-    #expand fit options
-    elastic, background, width, res_norm = Fit
-
-    #convert true/false to 1/0 for fortran
-    o_el = 1 if elastic else 0
-    o_w1 = 1 if width else 0
-    o_res = 1 if res_norm else 0
-
-    #fortran code uses background choices defined using the following numbers
-    if background == 'Sloping':
-        o_bgd = 2
-    elif background == 'Flat':
-        o_bgd = 1
-    elif background == 'Zero':
-        o_bgd = 0
-
-    fitOp = [o_el, o_bgd, o_w1, o_res]
-
-    workdir = config['defaultsave.directory']
-    if not os.path.isdir(workdir):
-        workdir = os.getcwd()
-        logger.information('Default Save directory is not set. Defaulting to current working Directory: ' + workdir)
-
-    array_len = 4096                           # length of array in Fortran
-    CheckXrange(erange,'Energy')
-    nbin,nrbin = nbins[0],nbins[1]
-    logger.information('Sample is ' + samWS)
-    logger.information('Resolution is ' + resWS)
-    CheckAnalysers(samWS,resWS)
-    nsam,ntc = CheckHistZero(samWS)
-
-    if Loop != True:
-        nsam = 1
-
-    efix = getEfixed(samWS)
-    theta,Q = GetThetaQ(samWS)
-    nres = CheckHistZero(resWS)[0]
-    if nres == 1:
-        prog = 'Qst'                        # res file
-    else:
-        raise ValueError('Stretched Exp ONLY works with RES file')
-    logger.information(' Number of spectra = '+str(nsam))
-    logger.information(' Erange : '+str(erange[0])+' to '+str(erange[1]))
-
-    fname = samWS[:-4] + '_'+ prog
-    wrks=os.path.join(workdir, samWS[:-4])
-    logger.information(' lptfile : ' + wrks +'_Qst.lpt')
-    lwrk=len(wrks)
-    wrks.ljust(140,' ')
-    wrkr=resWS
-    wrkr.ljust(140,' ')
-    Nbet,Nsig = nbs[0], nbs[1]
-    eBet0 = np.zeros(Nbet)                  # set errors to zero
-    eSig0 = np.zeros(Nsig)                  # set errors to zero
-    rscl = 1.0
-    Qaxis = ''
-    for m in range(0,nsam):
-        logger.information('Group ' +str(m)+ ' at angle '+ str(theta[m]))
-        nsp = m+1
-        nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(samWS,m,erange,nbin)
-        Ndat = nout[0]
-        Imin = nout[1]
-        Imax = nout[2]
-        Nb,Xb,Yb,_ = GetXYE(resWS,0,array_len)
-        numb = [nsam, nsp, ntc, Ndat, nbin, Imin, Imax, Nb, nrbin, Nbet, Nsig]
-        reals = [efix, theta[m], rscl, bnorm]
-        xsout,ysout,xbout,ybout,zpout=Que.quest(numb,Xv,Yv,Ev,reals,fitOp,\
-                                            Xdat,Xb,Yb,wrks,wrkr,lwrk)
-        dataXs = xsout[:Nsig]               # reduce from fixed Fortran array
-        dataYs = ysout[:Nsig]
-        dataXb = xbout[:Nbet]
-        dataYb = ybout[:Nbet]
-        zpWS = fname + '_Zp' +str(m)
-        if m > 0:
-            Qaxis += ','
-        Qaxis += str(Q[m])
-
-        dataXz = []
-        dataYz = []
-        dataEz = []
-
-        for n in range(0,Nsig):
-            yfit_list = np.split(zpout[:Nsig*Nbet],Nsig)
-            dataYzp = yfit_list[n]
-
-            dataXz = np.append(dataXz,xbout[:Nbet])
-            dataYz = np.append(dataYz,dataYzp[:Nbet])
-            dataEz = np.append(dataEz,eBet0)
-
-        CreateWorkspace(OutputWorkspace=zpWS, DataX=dataXz, DataY=dataYz, DataE=dataEz,
-                        Nspec=Nsig, UnitX='MomentumTransfer',
-                        VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=dataXs)
-
-        unitx = mtd[zpWS].getAxis(0).setUnit("Label")
-        unitx.setLabel('beta' , '')
-        unity = mtd[zpWS].getAxis(1).setUnit("Label")
-        unity.setLabel('sigma' , '')
-
-        if m == 0:
-            xSig = dataXs
-            ySig = dataYs
-            eSig = eSig0
-            xBet = dataXb
-            yBet = dataYb
-            eBet = eBet0
-            groupZ = zpWS
-        else:
-            xSig = np.append(xSig,dataXs)
-            ySig = np.append(ySig,dataYs)
-            eSig = np.append(eSig,eSig0)
-            xBet = np.append(xBet,dataXb)
-            yBet = np.append(yBet,dataYb)
-            eBet = np.append(eBet,eBet0)
-            groupZ = groupZ +','+ zpWS
-
-    #create workspaces for sigma and beta
-    CreateWorkspace(OutputWorkspace=fname+'_Sigma', DataX=xSig, DataY=ySig, DataE=eSig,\
-        Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
-    unitx = mtd[fname+'_Sigma'].getAxis(0).setUnit("Label")
-    unitx.setLabel('sigma' , '')
-
-    CreateWorkspace(OutputWorkspace=fname+'_Beta', DataX=xBet, DataY=yBet, DataE=eBet,\
-        Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
-    unitx = mtd[fname+'_Beta'].getAxis(0).setUnit("Label")
-    unitx.setLabel('beta' , '')
-
-    group = fname + '_Sigma,'+ fname + '_Beta'
-
-    fit_workspace = fname+'_Fit'
-    contour_workspace = fname+'_Contour'
-    GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fit_workspace)
-    GroupWorkspaces(InputWorkspaces=groupZ,OutputWorkspace=contour_workspace)
-
-    #add sample logs to the output workspaces
-    CopyLogs(InputWorkspace=samWS, OutputWorkspace=fit_workspace)
-    QuestAddSampleLogs(fit_workspace, resWS, background, elastic, erange, nbin, Nsig, Nbet)
-    CopyLogs(InputWorkspace=samWS, OutputWorkspace=contour_workspace)
-    QuestAddSampleLogs(contour_workspace, resWS, background, elastic, erange, nbin, Nsig, Nbet)
-
-    if Save:
-        fpath = os.path.join(workdir,fit_workspace+'.nxs')
-        SaveNexusProcessed(InputWorkspace=fit_workspace, Filename=fpath)
-
-        cpath = os.path.join(workdir,contour_workspace+'.nxs')
-        SaveNexusProcessed(InputWorkspace=contour_workspace, Filename=cpath)
-
-        logger.information('Output file for Fit : ' + fpath)
-        logger.information('Output file for Contours : ' + cpath)
-
-    if Plot != 'None' and Loop == True:
-        QuestPlot(fname,Plot)
-    EndTime('Quest')
-
-def QuestAddSampleLogs(workspace, res_workspace, background, elastic_peak, e_range, sample_binning, sigma, beta):
-    energy_min, energy_max = e_range
-
-    AddSampleLog(Workspace=workspace, LogName="res_file",
-                 LogType="String", LogText=res_workspace)
-    AddSampleLog(Workspace=workspace, LogName="background",
-                 LogType="String", LogText=str(background))
-    AddSampleLog(Workspace=workspace, LogName="elastic_peak",
-                 LogType="String", LogText=str(elastic_peak))
-    AddSampleLog(Workspace=workspace, LogName="energy_min",
-                 LogType="Number", LogText=str(energy_min))
-    AddSampleLog(Workspace=workspace, LogName="energy_max",
-                 LogType="Number", LogText=str(energy_max))
-    AddSampleLog(Workspace=workspace, LogName="sample_binning",
-                 LogType="Number", LogText=str(sample_binning))
-    AddSampleLog(Workspace=workspace, LogName="sigma",
-                 LogType="Number", LogText=str(sigma))
-    AddSampleLog(Workspace=workspace, LogName="beta",
-                 LogType="Number", LogText=str(beta))
-
-
-def QuestPlot(inputWS,Plot):
-    if Plot == 'Sigma' or Plot == 'All':
-        MTD_PLOT.importMatrixWorkspace(inputWS+'_Sigma').plotGraph2D()
-    if Plot == 'Beta' or Plot == 'All':
-        MTD_PLOT.importMatrixWorkspace(inputWS+'_Beta').plotGraph2D()
 
 # ResNorm programs
 def ResNormRun(vname,rname,erange,nbin,Plot='None',Save=False):