Add new integration algorithm for single-crystal Bragg peaks

Framework/PythonInterface/mantid/api/src/Exports/IPeakFunction.cpp

@@ -26,6 +26,7 @@ void export_IPeakFunction() {
            "Calculate the values of the function for the given x values. The "
            "output should be stored in the out array")
       .def("fwhm", &IPeakFunction::fwhm, arg("self"), "Returns the fwhm of the peak function.")
+      .def("setFwhm", &IPeakFunction::setFwhm, (arg("self"), arg("new_fwhm")), "Sets FWHM of peak.")
       .def("intensity", &IPeakFunction::intensity, arg("self"), "Returns the integral intensity of the peak function.")
       .def("intensityError", &IPeakFunction::intensityError, arg("self"),
            "Returns the integral intensity error of the peak function due to uncertainties in uncorrelated fit "
@@ -34,5 +35,5 @@ void export_IPeakFunction() {
            "Changes the integral intensity of the peak function by setting its "
            "height.")
       .def("setHeight", &IPeakFunction::setHeight, arg("self"), "Sets height of the peak function.")
-      .def("setCentre", &IPeakFunction::setCentre, arg("self"), "Sets setCentre of the peak function.");
+      .def("setCentre", &IPeakFunction::setCentre, arg("self"), "Sets the centre of the peak function.");
 }

Framework/PythonInterface/plugins/CMakeLists.txt

@@ -90,6 +90,7 @@ set(PYTHON_PLUGINS
     algorithms/IndexSatellitePeaks.py
     algorithms/IndirectTransmission.py
     algorithms/InelasticEMUauReduction.py
+    algorithms/IntegratePeaks1DProfile.py
     algorithms/IntegratePeaksProfileFitting.py
     algorithms/IntegratePeaksSkew.py
     algorithms/IntegratePeaksShoeboxTOF.py

Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py

@@ -366,7 +366,7 @@ def PyExec(self):
                 ipos_predicted = [peak_data.irow, peak_data.icol, ix]
                 det_edges = peak_data.det_edges if not integrate_on_edge else None
 
-                intens, sigma, status, ipos, nrows, ncols, nbins = integrate_peak(
+                intens, sigma, i_over_sig, status, ipos, nrows, ncols, nbins = integrate_peak(
                     ws,
                     peaks,
                     ipk,
@@ -383,7 +383,6 @@ def PyExec(self):
                     weak_peak_threshold,
                     do_optimise_shoebox,
                 )
-
                 if status == PEAK_STATUS.WEAK and do_optimise_shoebox and weak_peak_strategy == "NearestStrongPeak":
                     # look for possible strong peaks at any TOF in the window (won't know if strong until all pks integrated)
                     ipks_near, _ = find_ipks_in_window(ws, peaks, ispecs, ipk)
@@ -392,12 +391,9 @@ def PyExec(self):
                 else:
                     if status == PEAK_STATUS.STRONG:
                         ipks_strong.append(ipk)
-                    if output_file:
-                        # save result for plotting
-                        i_over_sig = intens / sigma if sigma > 0 else 0.0
-                        results[ipk] = ShoeboxResult(
-                            ipk, peak, x, y, [nrows, ncols, nbins], ipos, [peak_data.irow, peak_data.icol, ix], i_over_sig, status
-                        )
+                    results[ipk] = ShoeboxResult(
+                        ipk, peak, x, y, [nrows, ncols, nbins], ipos, [peak_data.irow, peak_data.icol, ix], i_over_sig, status
+                    )  # use this to get strong peak shoebox dimensions even if no plotting
                 # scale summed intensity by bin width to get integrated area
                 intens = intens * bin_width
                 sigma = sigma * bin_width
@@ -437,14 +433,13 @@ def PyExec(self):
                 # integrate at previously found ipos
                 ipos = [*np.argwhere(ispecs == weak_pk.ispec)[0], np.argmin(abs(x - weak_pk.tof))]
                 det_edges = peak_data.det_edges if not integrate_on_edge else None
-                intens, sigma, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges)
+                intens, sigma, i_over_sig, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges)
                 # scale summed intensity by bin width to get integrated area
                 intens = intens * weak_pk.tof_bin_width
                 sigma = sigma * weak_pk.tof_bin_width
                 set_peak_intensity(peak, intens, sigma, do_lorz_cor)
                 if output_file:
                     # save result for plotting
-                    i_over_sig = intens / sigma if sigma > 0 else 0.0
                     peak_shape = [nrows, ncols, nbins]
                     ipos_predicted = [peak_data.irow, peak_data.icol, np.argmin(abs(x - pk_tof))]
                     results[ipk] = ShoeboxResult(
@@ -500,23 +495,23 @@ def integrate_peak(
     ws, peaks, ipk, kernel, nrows, ncols, nbins, x, y, esq, ispecs, ipos_predicted, det_edges, weak_peak_threshold, do_optimise_shoebox
 ):
     # perform initial integration
-    intens_over_sig = convolve_shoebox(y, esq, kernel)
+    intens_over_sig = convolve_shoebox(y, esq, kernel)  # array of I/sigma same size as data
 
     # identify best shoebox position near peak
     ipos = find_nearest_peak_in_data_window(intens_over_sig, ispecs, x, ws, peaks, ipk, *ipos_predicted)
 
     # perform final integration if required
-    intens, sigma = 0.0, 0.0
+    intens, sigma, i_over_sig = 0.0, 0.0, 0.0
     status = PEAK_STATUS.NO_PEAK
     if ipos is not None:
         # integrate at that position (without smoothing I/sigma)
-        intens, sigma, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges)
+        intens, sigma, i_over_sig, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges)
         if status == PEAK_STATUS.STRONG and do_optimise_shoebox:
             ipos, (nrows, ncols, nbins) = optimise_shoebox(y, esq, (nrows, ncols, nbins), ipos)
             kernel = make_kernel(nrows, ncols, nbins)
             # re-integrate but this time check for overlap with edge
-            intens, sigma, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges)
-    return intens, sigma, status, ipos, nrows, ncols, nbins
+            intens, sigma, i_over_sig, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges)
+    return intens, sigma, i_over_sig, status, ipos, nrows, ncols, nbins
 
 
 def round_up_to_odd_number(number):
@@ -647,9 +642,9 @@ def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edge
     )
 
     status = PEAK_STATUS.NO_PEAK
+    intens, sigma, i_over_sig = 0.0, 0.0, 0.0
     if det_edges is not None and det_edges[slices[:-1]].any():
         status = PEAK_STATUS.ON_EDGE
-        intens, sigma = 0.0, 0.0
     else:
         if y[slices].size != kernel.size:
             # peak is partially on edge, but we continue to integrate with a partial kernel
@@ -661,7 +656,7 @@ def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edge
                 # get index in y where kernel starts
                 iy_start = ipos[idim] - kernel.shape[idim] // 2
                 if iy_start < 0:
-                    istart = -iy_start  # chop of ii elements at the begninning of the kernel along this dimension
+                    istart = -iy_start  # chop of ii elements at the beginning of the kernel along this dimension
                 elif iy_start > y.shape[idim] - kernel.shape[idim]:
                     iend = y.shape[idim] - iy_start  # include only up to number of elements remaining in y
                 kernel_slices.append(slice(istart, iend))
@@ -671,9 +666,10 @@ def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edge
             kernel[inegative] = -np.sum(kernel[~inegative]) / np.sum(inegative)
         intens = np.sum(y[slices] * kernel)
         sigma = np.sqrt(np.sum(esq[slices] * (kernel**2)))
+        i_over_sig = intens / sigma if sigma > 0 else 0.0
         if status == PEAK_STATUS.NO_PEAK:
-            status = PEAK_STATUS.STRONG if intens / sigma > weak_peak_threshold else PEAK_STATUS.WEAK
-    return intens, sigma, status
+            status = PEAK_STATUS.STRONG if i_over_sig > weak_peak_threshold else PEAK_STATUS.WEAK
+    return intens, sigma, i_over_sig, status
 
 
 def optimise_shoebox(y, esq, peak_shape, ipos, nfail_max=2):

Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt

@@ -69,6 +69,7 @@ set(TEST_PY_FILES
     IndexPeaksTest.py
     IndirectTransmissionTest.py
     IndexSatellitePeaksTest.py
+    IntegratePeaks1DProfileTest.py
     IntegratePeaksShoeboxTOFTest.py
     IntegratePeaksSkewTest.py
     LeadPressureCalcTest.py

Framework/PythonInterface/test/python/plugins/algorithms/IntegratePeaks1DProfileTest.py

@@ -0,0 +1,137 @@
+# Mantid Repository : https://github.com/mantidproject/mantid
+#
+# Copyright © 2023 ISIS Rutherford Appleton Laboratory UKRI,
+#   NScD Oak Ridge National Laboratory, European Spallation Source,
+#   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
+# SPDX - License - Identifier: GPL - 3.0 +
+import unittest
+from unittest import mock
+from mantid.simpleapi import (
+    IntegratePeaks1DProfile,
+    CreatePeaksWorkspace,
+    AddPeak,
+    Rebin,
+    LoadEmptyInstrument,
+    AnalysisDataService,
+    SortPeaksWorkspace,
+    CloneWorkspace,
+    BackToBackExponential,
+    DeleteTableRows,
+)
+import numpy as np
+import tempfile
+import shutil
+from os import path
+
+
+class IntegratePeaks1DProfileTest(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        # load empty instrument with RectangularDetector banks and create a peak table
+        cls.ws = LoadEmptyInstrument(InstrumentName="SXD", OutputWorkspace="SXD")
+        axis = cls.ws.getAxis(0)
+        axis.setUnit("TOF")
+        # rebin to get 20 TOF bins
+        cls.ws = Rebin(InputWorkspace=cls.ws, OutputWorkspace=cls.ws.name(), Params="10000,25,10500", PreserveEvents=False)
+        cls.ws += 50  # add constant background
+
+        cls.peaks_edge = CreatePeaksWorkspace(InstrumentWorkspace=cls.ws, NumberOfPeaks=0, OutputWorkspace="peaks_incl_edge")
+        cls.detids = [6049, 4220]
+        ispecs = cls.ws.getIndicesFromDetectorIDs(cls.detids)
+        # simulate peak
+        cls.pk_tof = 10250
+        pk_func = BackToBackExponential(I=1e4, A=0.9, X0=cls.pk_tof)  # override default A so cost-func not 0 in fit
+        for ipk, detid in enumerate(cls.detids):
+            AddPeak(PeaksWorkspace=cls.peaks_edge, RunWorkspace=cls.ws, TOF=cls.pk_tof, DetectorID=detid)
+            pk_func.function.setMatrixWorkspace(cls.ws, ispecs[ipk], 0.0, 0.0)
+            for ispec in np.arange(ispecs[ipk], ispecs[ipk] + 2):
+                y = cls.ws.dataY(int(ispec))
+                y += pk_func(cls.ws.readX(int(ispec))[:-1])  # note shifts peak centre by half a bin
+                cls.ws.setE(int(ispec), np.sqrt(y))
+
+        # clone workspace and delete edge peak so only 1 peak (quicker for most tests)
+        cls.peaks = CloneWorkspace(cls.peaks_edge, OutputWorkspace="peaks")
+        DeleteTableRows(cls.peaks, Rows=[1])
+
+        # default IntegratePeaks1DProfile kwargs
+        cls.profile_kwargs = {
+            "GetNBinsFromBackToBackParams": True,
+            "NFWHM": 6,
+            "CostFunction": "RSq",
+            "FixPeakParameters": "A",
+            "FractionalChangeDSpacing": 0.025,
+            "IntegrateIfOnEdge": True,
+            "LorentzCorrection": False,
+        }
+
+        # output file dir
+        cls._test_dir = tempfile.mkdtemp()
+
+    @classmethod
+    def tearDownClass(cls):
+        AnalysisDataService.clear()
+        shutil.rmtree(cls._test_dir)
+
+    def test_exec_IntegrateIfOnEdge_True(self):
+        out = IntegratePeaks1DProfile(
+            InputWorkspace=self.ws, PeaksWorkspace=self.peaks_edge, OutputWorkspace="peaks_int_1", **self.profile_kwargs
+        )
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 19.59, delta=1e-1)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[1], 19.59, delta=1e-1)
+
+    def test_exec_IntegrateIfOnEdge_False(self):
+        kwargs = self.profile_kwargs.copy()
+        kwargs["IntegrateIfOnEdge"] = False
+        out = IntegratePeaks1DProfile(InputWorkspace=self.ws, PeaksWorkspace=self.peaks_edge, OutputWorkspace="peaks_int_2", **kwargs)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 19.59, delta=1e-1)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[1], 0.0, delta=1e-2)
+
+    def test_exec_poisson_cost_func(self):
+        kwargs = self.profile_kwargs.copy()
+        kwargs["CostFunction"] = "Poisson"
+        out = IntegratePeaks1DProfile(InputWorkspace=self.ws, PeaksWorkspace=self.peaks, OutputWorkspace="peaks_int_3", **kwargs)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 19.60, delta=1e-2)
+
+    def test_exec_chisq_cost_func(self):
+        kwargs = self.profile_kwargs.copy()
+        kwargs["CostFunction"] = "ChiSq"
+        out = IntegratePeaks1DProfile(InputWorkspace=self.ws, PeaksWorkspace=self.peaks, OutputWorkspace="peaks_int_4", **kwargs)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 19.60, delta=1e-2)
+
+    def test_exec_hessian_error_strategy(self):
+        kwargs = self.profile_kwargs.copy()
+        kwargs["ErrorStrategy"] = "Hessian"
+        out = IntegratePeaks1DProfile(InputWorkspace=self.ws, PeaksWorkspace=self.peaks, OutputWorkspace="peaks_int_5", **kwargs)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 57922.22, delta=1e-2)  # not realistic fit
+
+    def test_exec_IOverSigmaThreshold_respected(self):
+        kwargs = self.profile_kwargs.copy()
+        kwargs["IOverSigmaThreshold"] = 100  # set this higher than I/sigma in any pixel i.e. should not fit any pixels
+        out = IntegratePeaks1DProfile(InputWorkspace=self.ws, PeaksWorkspace=self.peaks, OutputWorkspace="peaks_int_6", **kwargs)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 0.0, delta=1e-2)
+
+    def test_exec_gaussian_peak_func(self):
+        kwargs = self.profile_kwargs.copy()
+        kwargs["PeakFunction"] = "Gaussian"
+        kwargs["FixPeakParameters"] = ""
+        out = IntegratePeaks1DProfile(InputWorkspace=self.ws, PeaksWorkspace=self.peaks, OutputWorkspace="peaks_int_7", **kwargs)
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 19.64, delta=1e-2)
+
+    def test_exec_OutputFile(self):
+        out_file = path.join(self._test_dir, "out.pdf")
+        IntegratePeaks1DProfile(
+            InputWorkspace=self.ws, PeaksWorkspace=self.peaks, OutputWorkspace="peaks_int_8", OutputFile=out_file, **self.profile_kwargs
+        )
+        # check output file saved
+        self.assertTrue(path.exists(out_file))
+
+    def test_exec_FractionalChangeDSpacing(self):
+        kwargs = self.profile_kwargs.copy()
+        kwargs["FractionalChangeDSpacing"] = 1e-8
+        out = IntegratePeaks1DProfile(InputWorkspace=self.ws, PeaksWorkspace=self.peaks, OutputWorkspace="peaks_int_9", **kwargs)
+        # I/sigma different as center constrained
+        self.assertAlmostEqual(out.column("Intens/SigInt")[0], 20.41, delta=1e-2)
+
+
+if __name__ == "__main__":
+    unittest.main()

Testing/SystemTests/tests/framework/SXDAnalysis.py

@@ -135,38 +135,67 @@ def validate(self):
         self.tolerance_is_rel_err = True
         return self.integrated_peaks, "SXD23767_found_peaks_integrated.nxs"
 
-    class SXDIntegrateDataShoebox(systemtesting.MantidSystemTest):
-        def cleanup(self):
-            ADS.clear()
-
-        def runTest(self):
-            sxd = SXD(vanadium_runno=23769, empty_runno=23768)
-            # load data and convert to Qlab
-            ws = LoadNexus(Filename="SXD23767_processed.nxs", OutputWorkspace="SXD23767_processed")
-            runno = 23767
-            sxd.set_ws(runno, ws)
-            # create peak table
-            peaks = CreatePeaksWorkspace(InstrumentWorkspace=ws, NumberOfPeaks=0, OutputWorkspace="peaks")
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=8303.3735339704781, DetectorID=7646)
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=19166.422281671705, DetectorID=20009)
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=1582.4728457415549, DetectorID=25111)
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=1734.9661945986509, DetectorID=28271)
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=5241.7309709929505, DetectorID=1380)
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=1665.2648421691604, DetectorID=8176)
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=7690.5100117616385, DetectorID=22799)
-            AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=8289.2655006925579, DetectorID=7327)
-            sxd.set_peaks(runno, peaks, PEAK_TYPE.FOUND)
-            # integrate
-            sxd.integrate_data(
-                INTEGRATION_TYPE.SHOEBOX,
-                PEAK_TYPE.FOUND,
-                GetNBinsFromBackToBackParams=True,
-                WeakPeakThreshold=20.0,
-                LorentzCorrection=False,
-                WeakPeakStrategy="NearestStrongPeak",
-            )
-            self.integrated_peaks = sxd.get_peaks_name(runno, PEAK_TYPE.FOUND, INTEGRATION_TYPE.SHOEBOX)
-
-        def validate(self):
-            intens_over_sigma = [98.97, 0.0, 10.5, 10.87, 135.83, 0.0, -0.4, 1.08]
-            self.assertTrue(np.allclose(np.round(self.integrated_peaks.column("Intens/SigInt"), 2), intens_over_sigma))
+
+class SXDIntegrateDataShoebox(systemtesting.MantidSystemTest):
+    def cleanup(self):
+        ADS.clear()
+
+    def runTest(self):
+        sxd = SXD(vanadium_runno=23769, empty_runno=23768)
+        runno = 23767
+        _load_data_and_add_peaks(sxd, runno)
+
+        # integrate
+        sxd.integrate_data(
+            INTEGRATION_TYPE.SHOEBOX,
+            PEAK_TYPE.FOUND,
+            GetNBinsFromBackToBackParams=True,
+            WeakPeakThreshold=20.0,
+            LorentzCorrection=False,
+            WeakPeakStrategy="NearestStrongPeak",
+        )
+        self.integrated_peaks = sxd.get_peaks(runno, PEAK_TYPE.FOUND, INTEGRATION_TYPE.SHOEBOX)
+
+    def validate(self):
+        intens_over_sigma = [0.0, 10.50, 9.56, 135.83, 0.0]
+        self.assertTrue(np.allclose(self.integrated_peaks.column("Intens/SigInt"), intens_over_sigma, atol=1e-2))
+
+
+class SXDIntegrateData1DProfile(systemtesting.MantidSystemTest):
+    def cleanup(self):
+        ADS.clear()
+
+    def runTest(self):
+        sxd = SXD(vanadium_runno=23769, empty_runno=23768)
+        runno = 23767
+        _load_data_and_add_peaks(sxd, runno)
+
+        # integrate
+        sxd.integrate_data(
+            INTEGRATION_TYPE.PROFILE,
+            PEAK_TYPE.FOUND,
+            GetNBinsFromBackToBackParams=True,
+            NFWHM=6,
+            CostFunction="RSq",
+            FixPeakParameters="A",
+            FractionalChangeDSpacing=0.025,
+            IntegrateIfOnEdge=False,
+        )
+        self.integrated_peaks = sxd.get_peaks(runno, PEAK_TYPE.FOUND, INTEGRATION_TYPE.PROFILE)
+
+    def validate(self):
+        intens_over_sigma = [0.0, 18.45, 17.11, 134.29, 0.0]
+        self.assertTrue(np.allclose(self.integrated_peaks.column("Intens/SigInt"), intens_over_sigma, atol=1e-2))
+
+
+def _load_data_and_add_peaks(sxd, runno):
+    ws = LoadNexus(Filename="SXD23767_processed.nxs", OutputWorkspace="SXD23767_processed")
+    sxd.set_ws(runno, ws)
+    # create peak table
+    peaks = CreatePeaksWorkspace(InstrumentWorkspace=ws, NumberOfPeaks=0, OutputWorkspace="peaks")
+    AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=19166.422281671705, DetectorID=20009)
+    AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=1582.4728457415549, DetectorID=25111)
+    AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=1734.9661945986509, DetectorID=28271)
+    AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=5241.7309709929505, DetectorID=1380)
+    AddPeak(PeaksWorkspace="peaks", RunWorkspace=ws, TOF=1665.2648421691604, DetectorID=8176)
+    sxd.set_peaks(runno, peaks, PEAK_TYPE.FOUND)