Fix FitIncidentSpectrum fitting wrong data by loading raw.

Also fixes FitIncidentSpectrum.py and CalculatePlaczekSelfScattering so
that output workspace has same x units

Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp

@@ -5,6 +5,7 @@
 //     & Institut Laue - Langevin
 // SPDX - License - Identifier: GPL - 3.0 +
 #include "MantidAlgorithms/CalculatePlaczekSelfScattering.h"
+#include "MantidAPI/Axis.h"
 #include "MantidAPI/Sample.h"
 #include "MantidAPI/WorkspaceFactory.h"
 #include "MantidDataObjects/Workspace2D.h"
@@ -12,6 +13,7 @@
 #include "MantidGeometry/Instrument/DetectorInfo.h"
 #include "MantidKernel/Atom.h"
 #include "MantidKernel/Material.h"
+#include "MantidKernel/Unit.h"
 
 #include <utility>
 
@@ -164,11 +166,12 @@ void CalculatePlaczekSelfScattering::exec() {
       }
     }
   }
+  std::string unit = inWS->getAxis(0)->unit()->unitID();
   Mantid::API::Algorithm_sptr childAlg =
       createChildAlgorithm("CreateWorkspace");
   childAlg->setProperty("DataX", xLambdas);
   childAlg->setProperty("DataY", placzekCorrection);
-  childAlg->setProperty("UnitX", "Wavelength");
+  childAlg->setProperty("UnitX", unit);
   childAlg->setProperty("NSpec", nSpec);
   childAlg->setProperty("ParentWorkspace", inWS);
   childAlg->setProperty("Distribution", true);

Framework/PythonInterface/plugins/algorithms/FitIncidentSpectrum.py

@@ -44,6 +44,11 @@ class FitIncidentSpectrum(PythonAlgorithm):
                                     direction=Direction.Output),
             doc='Output workspace containing the fit and it\'s first derivative.')
 
+        self.declareProperty(
+            name='SpectrumIndex',
+            defaultValue=0,
+            doc='Workspace index of the spectra to be fitted (Defaults to the first index.)')
+
         self.declareProperty(FloatArrayProperty(name="BinningForCalc",
                                                 validator=RebinParamsValidator(AllowEmpty=True),
                                                 direction=Direction.Input),
@@ -67,29 +72,28 @@ class FitIncidentSpectrum(PythonAlgorithm):
     def _setup(self):
         self._input_ws = self.getProperty('InputWorkspace').value
         self._output_ws = self.getProperty('OutputWorkspace').valueAsStr
+        self._incident_index = self.getProperty('SpectrumIndex').value
+        self._binning_for_calc = self.getProperty('BinningForCalc').value
         self._binning_for_fit = self.getProperty('BinningForFit').value
         self._fit_spectrum_with = self.getProperty('FitSpectrumWith').value
-        self._binning_for_calc = self.getProperty('BinningForCalc').value
-        if not self._binning_for_calc.all():
+        if self._binning_for_calc.size == 0:
             x = self._input_ws.readX(0)
-            self._binning_for_calc = [str(i) for i in [min(x), x[1] - x[0], max(x) + x[1] - x[0]]]
+            self._binning_for_calc = [i for i in [min(x), x[1] - x[0], max(x) + x[1] - x[0]]]
 
     def PyExec(self):
         self._setup()
-
         x = np.arange(self._binning_for_calc[0], self._binning_for_calc[2], self._binning_for_calc[1])
-        incident_index = 0
         if self._binning_for_fit.size == 0:
-            x_fit = np.array(self._input_ws.readX(incident_index))
-            y_fit = np.array(self._input_ws.readY(incident_index))
+            x_fit = np.array(self._input_ws.readX(self._incident_index))
+            y_fit = np.array(self._input_ws.readY(self._incident_index))
         else:
             rebinned = Rebin(
                 self._input_ws,
                 Params=self._binning_for_fit,
                 PreserveEvents=True,
                 StoreInADS=False)
-            x_fit = np.array(rebinned.readX(incident_index))
-            y_fit = np.array(rebinned.readY(incident_index))
+            x_fit = np.array(rebinned.readX(self._incident_index))
+            y_fit = np.array(rebinned.readY(self._incident_index))
 
         if len(x_fit) != len(y_fit):
             x_fit = x_fit[:-1]
@@ -110,10 +114,11 @@ class FitIncidentSpectrum(PythonAlgorithm):
             fit, fit_prime = self.fit_cubic_spline_with_gauss_conv(x_fit, y_fit, x, sigma=0.5)
 
         # Create output workspace
+        unit = self._input_ws.getAxis(0).getUnit().unitID()
         output_workspace = CreateWorkspace(
             DataX=x,
             DataY=np.append(fit, fit_prime),
-            UnitX='Wavelength',
+            UnitX=unit,
             NSpec=2,
             Distribution=False,
             ParentWorkspace=self._input_ws,

scripts/Diffraction/isis_powder/polaris.py

@@ -58,7 +58,8 @@ class Polaris(AbstractInst):
                                                   focus_file_path=focus_file_path,
                                                   merge_banks=self._inst_settings.merge_banks,
                                                   q_lims=q_lims,
-                                                  cal_file_name=cal_file_name)
+                                                  cal_file_name=cal_file_name,
+                                                  sample_details=self._sample_details)
         return pdf_output
 
     def set_sample_details(self, **kwargs):

scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py

@@ -82,13 +82,13 @@ def save_unsplined_vanadium(vanadium_ws, output_path):
         mantid.SaveNexus(InputWorkspace=vanadium_ws, Filename=output_path, Append=False)
 
 
-def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None, cal_file_name=None):
+def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None, cal_file_name=None, sample_details=None):
     focused_ws = _obtain_focused_run(run_number, focus_file_path)
-    focused_ws = mantid.ConvertUnits(InputWorkspace=focused_ws.name(), Target="MomentumTransfer")
 
     if merge_banks:
-        pdf_output = _generate_grouped_ts_pdf(focused_ws, q_lims, cal_file_name)
+        pdf_output = _generate_grouped_ts_pdf(run_number, focused_ws, q_lims, cal_file_name, sample_details)
     else:
+        focused_ws = mantid.ConvertUnits(InputWorkspace=focused_ws.name(), Target="MomentumTransfer")
         pdf_output = mantid.PDFFourierTransform(Inputworkspace=focused_ws, InputSofQType="S(Q)", PDFType="G(r)",
                                                 Filter=True)
         pdf_output = mantid.RebinToWorkspace(WorkspaceToRebin=pdf_output, WorkspaceToMatch=pdf_output[4],
@@ -123,58 +123,42 @@ def _obtain_focused_run(run_number, focus_file_path):
     return focused_ws
 
 
-def _generate_grouped_ts_pdf(focused_ws, q_lims, cal_file_name):
-    group_bin_min = None
-    group_bin_max = None
-    group_bin_width = None
+def _generate_grouped_ts_pdf(run_number, focused_ws, q_lims, cal_file_name, sample_details):
+    raw_ws = mantid.LoadRaw(Filename='POL'+str(run_number))
+    mantid.SetSample(InputWorkspace=raw_ws,
+                     Geometry=common.generate_sample_geometry(sample_details),
+                     Material=common.generate_sample_material(sample_details))
+    fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=raw_ws,
+                                             SpectrumIndex=11,
+                                             FitSpectrumWith="GaussConvCubicSpline")
+    placzek_self_scattering = mantid.CalculatePlaczekSelfScattering(InputWorkspace=fit_spectra)
+    cal_workspace = mantid.LoadCalFile(InputWorkspace=placzek_self_scattering,
+                                       CalFileName=cal_file_name,
+                                       Workspacename='cal_workspace',
+                                       MakeOffsetsWorkspace=False,
+                                       MakeMaskWorkspace=False)
     for i in range(focused_ws.getNumberOfEntries()):
-        x_array = focused_ws.getItem(i).readX(0)
-        bin_min = x_array[0]
-        bin_max = x_array[-1]
-        bin_width = (x_array[-1] - x_array[0]) / x_array.size
-        binning = [bin_min, bin_width, bin_max]
-        if not group_bin_min:
-            group_bin_min = bin_min
-            group_bin_max = bin_max
-            group_bin_width = bin_width
-        else:
-            group_bin_min = min(group_bin_min, bin_min)
-            group_bin_max = max(group_bin_max, bin_max)
-            group_bin_width = min(group_bin_width, bin_width)
-        fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=focused_ws.getItem(i),
-                                                 BinningForFit=binning,
-                                                 BinningForCalc=binning,
-                                                 FitSpectrumWith="GaussConvCubicSpline")
-        # find the flattest part of the fit spectra for background subtraction
-        lambda_prime = fit_spectra.dataY(1)
-        n_bg = int(lambda_prime.size*0.05)
-        idx = np.argpartition(np.absolute(lambda_prime), n_bg)
-        background = np.mean(fit_spectra.dataY(0)[idx[:n_bg]])
-        placzek_self_scattering = mantid.CalculatePlaczekSelfScattering(InputWorkspace=fit_spectra)
-        cal_workspace = mantid.LoadCalFile(InputWorkspace=placzek_self_scattering,
-                                           CalFileName=cal_file_name,
-                                           Workspacename='cal_workspace',
-                                           MakeOffsetsWorkspace=False,
-                                           MakeMaskWorkspace=False)
         focused_correction = None
+        axis = focused_ws.getItem(i).getAxis(0)
+        unit = axis.getUnit().unitID()
         for spec_index in range(placzek_self_scattering.getNumberHistograms()):
             if cal_workspace.dataY(spec_index)[0] == i + 1:
                 if focused_correction is None:
-                    focused_correction = placzek_self_scattering.dataY(spec_index) + background
-                else:
-                    focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
+                    focused_correction = placzek_self_scattering.dataY(spec_index)
+                #else:
+                    #focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
         focused_correction_ws = mantid.CreateWorkspace(DataX=placzek_self_scattering.dataX(0),
                                                        DataY=focused_correction,
                                                        Distribution=True,
-                                                       UnitX='MomentumTransfer')
-        mantid.Rebin(InputWorkspace=focused_ws.getItem(i), OutputWorkspace=focused_ws.getItem(i), Params=binning)
-        focused_correction_ws = mantid.Rebin(InputWorkspace=focused_correction_ws, Params=binning)
+                                                       UnitX=unit)
+        mantid.RebinToWorkspace(WorkspaceToRebin=focused_correction_ws,
+                                WorkspaceToMatch=focused_ws.getItem(i),
+                                OutputWorkspace=focused_correction_ws)
         mantid.Subtract(LHSWorkspace=focused_ws.getItem(i),
                         RHSWorkspace=focused_correction_ws,
                         OutputWorkspace=focused_ws.getItem(i))
-    binning = [group_bin_min, group_bin_width, group_bin_max]
-    focused_data = mantid.Rebin(InputWorkspace=focused_ws, Params=binning)
-    focused_data_combined = mantid.ConjoinSpectra(InputWorkspaces=focused_data)
+    focused_data_combined = mantid.ConjoinSpectra(InputWorkspaces=focused_ws)
+    focused_data_combined = mantid.ConvertUnits(InputWorkspace=focused_data_combined, Target="MomentumTransfer")
     mantid.MatchSpectra(InputWorkspace=focused_data_combined,
                         OutputWorkspace=focused_data_combined,
                         ReferenceSpectrum=5)
@@ -195,27 +179,28 @@ def _generate_grouped_ts_pdf(focused_ws, q_lims, cal_file_name):
         q_max = q_lims[1, :]
     else:
         raise RuntimeError("q_lims is not valid")
-    pdf_x_array = focused_data_combined.readX(0)
+    bin_width = np.inf
     for i in range(q_min.size):
+        pdf_x_array = focused_data_combined.readX(i)
         q_min[i] = pdf_x_array[np.amin(np.where(pdf_x_array >= q_min[i]))]
         q_max[i] = pdf_x_array[np.amax(np.where(pdf_x_array <= q_max[i]))]
-    bin_width = pdf_x_array[1] - pdf_x_array[0]
+        bin_width = min(pdf_x_array[1] - pdf_x_array[0], bin_width)
     focused_data_combined = mantid.CropWorkspaceRagged(InputWorkspace=focused_data_combined, XMin=q_min, XMax=q_max)
     focused_data_combined = mantid.Rebin(InputWorkspace=focused_data_combined,
                                          Params=[min(q_min), bin_width, max(q_max)])
     focused_data_combined = mantid.SumSpectra(InputWorkspace=focused_data_combined,
                                               WeightedSum=True,
-                                              MultiplyBySpectra=False)
+                                             MultiplyBySpectra=False)
     pdf_output = mantid.PDFFourierTransform(Inputworkspace=focused_data_combined,
                                             InputSofQType="S(Q)",
                                             PDFType="G(r)",
                                             Filter=True)
 
-    common.remove_intermediate_workspace(cal_workspace)
-    common.remove_intermediate_workspace(fit_spectra)
-    common.remove_intermediate_workspace(focused_correction_ws)
-    common.remove_intermediate_workspace(focused_data)
-    common.remove_intermediate_workspace(placzek_self_scattering)
+    #common.remove_intermediate_workspace(cal_workspace)
+    #common.remove_intermediate_workspace(fit_spectra)
+    #common.remove_intermediate_workspace(focused_correction_ws)
+    #common.remove_intermediate_workspace(placzek_self_scattering)
+    #common.remove_intermediate_workspace(raw_ws)
     return pdf_output