Create basic workflow for generate_ts_pdf with merge_banks

Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp

@@ -32,7 +32,7 @@ getSampleSpeciesInfo(const API::MatrixWorkspace_const_sptr ws) {
   for (const Kernel::Material::Material::FormulaUnit element : formula) {
     const std::map<std::string, double> atomMap{
         {"mass", element.atom->mass},
-        {"stoich", element.multiplicity},
+        {"concentration", element.multiplicity},
         {"bSqrdBar", bSqrdBar}};
     atomSpecies[element.atom->symbol] = atomMap;
     totalStoich += element.multiplicity;
@@ -127,7 +127,9 @@ void CalculatePlaczekSelfScattering::exec() {
   const Geometry::DetectorInfo detInfo = inWS->detectorInfo();
   for (size_t detIndex = 0; detIndex < detInfo.size(); detIndex++) {
     if (!detInfo.isMonitor(detIndex)) {
+      if (!detInfo.l2(detIndex) == 0) {
       nReserve += 1;
+      }
     }
   }
   xLambdas.reserve(nReserve);
@@ -136,21 +138,30 @@ void CalculatePlaczekSelfScattering::exec() {
   int nSpec = 0;
   for (size_t detIndex = 0; detIndex < detInfo.size(); detIndex++) {
     if (!detInfo.isMonitor(detIndex)) {
-      const double pathLength = detInfo.l1() + detInfo.l2(detIndex);
-      const double f = detInfo.l1() / pathLength;
-      const double sinThetaBy2 = sin(detInfo.twoTheta(detIndex) / 2.0);
-      for (size_t xIndex = 0; xIndex < xLambda.size() - 1; xIndex++) {
-        const double term1 = (f - 1.0) * phi1[xIndex];
-        const double term2 = f * eps1[xIndex];
-        const double term3 = f - 3;
-        const double inelasticPlaczekSelfCorrection =
-            2.0 * (term1 - term2 + term3) * sinThetaBy2 * sinThetaBy2 *
-            summationTerm;
-        xLambdas.push_back(xLambda[xIndex]);
-        placzekCorrection.push_back(inelasticPlaczekSelfCorrection);
+      if (!detInfo.l2(detIndex) == 0) {
+        const double pathLength = detInfo.l1() + detInfo.l2(detIndex);
+        const double f = detInfo.l1() / pathLength;
+        const double sinThetaBy2 = sin(detInfo.twoTheta(detIndex) / 2.0);
+        for (size_t xIndex = 0; xIndex < xLambda.size() - 1; xIndex++) {
+          const double term1 = (f - 1.0) * phi1[xIndex];
+          const double term2 = f * eps1[xIndex];
+          const double term3 = f - 3;
+          const double inelasticPlaczekSelfCorrection =
+              2.0 * (term1 - term2 + term3) * sinThetaBy2 * sinThetaBy2 *
+              summationTerm;
+          xLambdas.push_back(xLambda[xIndex]);
+          placzekCorrection.push_back(inelasticPlaczekSelfCorrection);
+        }
+        xLambdas.push_back(xLambda.back());
+        nSpec += 1;
+      } else {
+        for (size_t xIndex = 0; xIndex < xLambda.size() - 1; xIndex++) {
+          xLambdas.push_back(xLambda[xIndex]);
+          placzekCorrection.push_back(0);
+        }
+        xLambdas.push_back(xLambda.back());
+        nSpec += 1;
       }
-      xLambdas.push_back(xLambda.back());
-      nSpec += 1;
     }
   }
   Mantid::API::Algorithm_sptr childAlg =

scripts/Diffraction/isis_powder/polaris.py

@@ -53,13 +53,12 @@ class Polaris(AbstractInst):
         # Generate pdf
         run_details = self._get_run_details(self._inst_settings.run_number)
         focus_file_path = self._generate_out_file_paths(run_details)["nxs_filename"]
-        run_details = self._get_run_details(run_number_string=self._inst_settings.run_number)
+        cal_file_name = self._inst_settings.calibration_dir + '/' + self._inst_settings.grouping_file_name
         pdf_output = polaris_algs.generate_ts_pdf(run_number=self._inst_settings.run_number,
                                                   focus_file_path=focus_file_path,
                                                   merge_banks=self._inst_settings.merge_banks,
-                                                  # q_lims=self._inst_settings.q_lims,
-                                                  # grouping_file_path=run_details.grouping_file_path
-                                                  )
+                                                  q_lims=self._inst_settings.q_lims,
+                                                  cal_file_name=cal_file_name)
         return pdf_output
 
     def set_sample_details(self, **kwargs):

scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py

@@ -5,7 +5,9 @@
 #     & Institut Laue - Langevin
 # SPDX - License - Identifier: GPL - 3.0 +
 from __future__ import (absolute_import, division, print_function)
+import numpy as np
 
+from mantid.api import WorkspaceGroup
 import mantid.simpleapi as mantid
 
 from isis_powder.routines import absorb_corrections, common
@@ -81,15 +83,58 @@ 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):
+def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None, cal_file_name=None):
     focused_ws = _obtain_focused_run(run_number, focus_file_path)
     pdf_output = mantid.ConvertUnits(InputWorkspace=focused_ws.name(), Target="MomentumTransfer")
 
     if merge_banks:
-        placzek_self_scattering = mantid.CalculatePlaczekSelfScattering(InputWorkspace=pdf_output)
-        pdf_output = mantid.Subtract(LHSWorkspace=pdf_output, RHSWorkspace=placzek_self_scattering)
-
-        pdf_output = mantid.MatchSpectra(InputWorkspace=pdf_output, ReferenceSpectrum=1)
+        group_bin_min = None
+        group_bin_max = None
+        group_bin_width = None
+        for i in range(pdf_output.getNumberOfEntries()):
+            x_array = pdf_output.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 = min(group_bin_max, bin_max)
+                group_bin_width = min(group_bin_width, bin_width)
+            fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=pdf_output.getItem(i),
+                                                     BinningForFit=binning,
+                                                     BinningForCalc=binning,
+                                                     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)
+            focused_correction = None
+            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)
+                    else:
+                        focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
+            mantid.CreateWorkspace(DataX=placzek_self_scattering.dataX(0),
+                                   DataY=focused_correction,
+                                   Distribution=True,
+                                   UnitX='MomentumTransfer',
+                                   OutputWorkspace='focused_correction_ws')
+            mantid.Rebin(InputWorkspace=pdf_output.getItem(i), OutputWorkspace='pdf_rebined', Params=binning)
+            mantid.Rebin(InputWorkspace='focused_correction_ws', OutputWorkspace='focused_correction_ws', Params=binning)
+            mantid.Subtract(LHSWorkspace='pdf_rebined', RHSWorkspace='focused_correction_ws', OutputWorkspace=pdf_output.getItem(i))
+        binning = [group_bin_min, group_bin_width, group_bin_max]
+        pdf_output = mantid.Rebin(InputWorkspace=pdf_output, Params=binning)
+        pdf_output_combined = mantid.ConjoinSpectra(InputWorkspaces='pdf_output')
+        mantid.MatchSpectra(InputWorkspace=pdf_output_combined, OutputWorkspace=pdf_output_combined, ReferenceSpectrum=1)
         if type(q_lims) == str:
             q_min = []
             q_max = []
@@ -102,16 +147,22 @@ def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None)
                         q_max.append(value_list[3])
             except IOError:
                 raise RuntimeError("q_lims is not valid")
-        elif type(q_lims) == list:
+        elif type(q_lims) == list or type(q_lims) == np.ndarray:
             q_min = q_lims[0, :]
             q_max = q_lims[1, :]
         else:
             raise RuntimeError("q_lims is not valid")
-        pdf_x_array = pdf_output.readX()
+        pdf_x_array = pdf_output_combined.readX(0)
+        for i in range(q_min.size):
+            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]
-        pdf_output = mantid.CropWorkspaceRagged(InputWorkspace=pdf_output, XMin=q_min, XMax=q_max)
-        pdf_output = mantid.Rebin(InputWorkspace=pdf_output, Params=[q_min, bin_width, q_max])
-        pdf_output = mantid.SumSpectra(InputWorkspace=pdf_output, WeightedSum=True)
+        pdf_output_combined = mantid.CropWorkspaceRagged(InputWorkspace=pdf_output_combined, XMin=q_min, XMax=q_max)
+        pdf_output_combined = mantid.Rebin(InputWorkspace=pdf_output_combined, Params=[min(q_min), bin_width, max(q_max)])
+        pdf_output_combined = mantid.SumSpectra(InputWorkspace=pdf_output_combined, WeightedSum=True, MultiplyBySpectra=False)
+        pdf_fourier_transform = mantid.PDFFourierTransform(Inputworkspace=pdf_output_combined, InputSofQType="S(Q)",
+                                                           PDFType="G(r)", Filter=True)
+        return pdf_fourier_transform
 
     pdf_output = mantid.PDFFourierTransform(Inputworkspace=pdf_output, InputSofQType="S(Q)", PDFType="G(r)",
                                             Filter=True)
@@ -120,6 +171,86 @@ def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None)
     return pdf_output
 
 
+def _generate_grouped_ts_pdf(focused_data, q_lims, cal_file_name):
+    group_bin_min = None
+    group_bin_max = None
+    group_bin_width = None
+    for i in range(focused_data.getNumberOfEntries()):
+        x_array = focused_data.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 = min(group_bin_max, bin_max)
+            group_bin_width = min(group_bin_width, bin_width)
+        fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=focused_data.getItem(i),
+                                                 BinningForFit=binning,
+                                                 BinningForCalc=binning,
+                                                 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)
+        focused_correction = None
+        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)
+                else:
+                    focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
+        mantid.CreateWorkspace(DataX=placzek_self_scattering.dataX(0),
+                               DataY=focused_correction,
+                               Distribution=True,
+                               UnitX='MomentumTransfer',
+                               OutputWorkspace='focused_correction_ws')
+        mantid.Rebin(InputWorkspace=focused_data.getItem(i), OutputWorkspace='focused_rebined', Params=binning)
+        mantid.Rebin(InputWorkspace='focused_correction_ws', OutputWorkspace='focused_correction_ws', Params=binning)
+        mantid.Subtract(LHSWorkspace='focused_rebined', RHSWorkspace='focused_correction_ws',
+                        OutputWorkspace=focused_data.getItem(i))
+    binning = [group_bin_min, group_bin_width, group_bin_max]
+    focused_data = mantid.Rebin(InputWorkspace=focused_data, Params=binning)
+    pdf_output_combined = mantid.ConjoinSpectra(InputWorkspaces=focused_data)
+    mantid.MatchSpectra(InputWorkspace=pdf_output_combined, OutputWorkspace=pdf_output_combined, ReferenceSpectrum=1)
+    if type(q_lims) == str:
+        q_min = []
+        q_max = []
+        try:
+            with open(q_lims, 'r') as f:
+                line_list = [line.rstrip('\n') for line in f]
+                for line in line_list[:-1]:
+                    value_list = line.split()
+                    q_min.append(value_list[2])
+                    q_max.append(value_list[3])
+        except IOError:
+            raise RuntimeError("q_lims is not valid")
+    elif type(q_lims) == list or type(q_lims) == np.ndarray:
+        q_min = q_lims[0, :]
+        q_max = q_lims[1, :]
+    else:
+        raise RuntimeError("q_lims is not valid")
+    pdf_x_array = pdf_output_combined.readX(0)
+    for i in range(q_min.size):
+        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]
+    focused_data_combined = mantid.CropWorkspaceRagged(InputWorkspace=pdf_output_combined, XMin=q_min, XMax=q_max)
+    focused_data_combined = mantid.Rebin(InputWorkspace=pdf_output_combined, Params=[min(q_min), bin_width, max(q_max)])
+    focused_data_combined = mantid.SumSpectra(InputWorkspace=pdf_output_combined, WeightedSum=True,
+                                            MultiplyBySpectra=False)
+    pdf_output = mantid.PDFFourierTransform(Inputworkspace=pdf_output_combined, InputSofQType="S(Q)",
+                                                       PDFType="G(r)", Filter=True)
+    return pdf_output
+
+
 def _obtain_focused_run(run_number, focus_file_path):
     """
     Searches for the focused workspace to use (based on user specified run number) in the ADS and then the output

scripts/Diffraction/isis_powder/polaris_routines/polaris_param_mapping.py

@@ -28,6 +28,7 @@ attr_mapping = [
     ParamMapEntry(ext_name="mode", int_name="mode", enum_class=POLARIS_CHOPPER_MODES,
                   optional=True),
     ParamMapEntry(ext_name="multiple_scattering", int_name="multiple_scattering", optional=True),
+    ParamMapEntry(ext_name="q_lims", int_name="q_lims"),
     ParamMapEntry(ext_name="raw_data_cropping_values", int_name="raw_data_crop_values"),
     ParamMapEntry(ext_name="run_number", int_name="run_number"),
     ParamMapEntry(ext_name="sample_empty", int_name="sample_empty", optional=True),