Add PeakShapeDetectorBin to IntegratePeaks1DProfile peak integration algorithm (#39040)

* DetectorBin peak shape added

* unit tests and release notes added

Author: warunawickramasingha

Framework/PythonInterface/plugins/algorithms/IntegratePeaks1DProfile.py

@@ -28,6 +28,7 @@
     PropertyCriterion,
     StringArrayProperty,
     UnitParams,
+    SpecialCoordinateSystem,
 )
 from mantid.dataobjects import Workspace2D
 from mantid.fitfunctions import FunctionWrapper, CompositeFunctionWrapper
@@ -42,6 +43,7 @@
 )
 from enum import Enum
 from typing import Sequence, Tuple
+from mantid.dataobjects import PeakShapeDetectorBin
 
 
 class PEAK_STATUS(Enum):
@@ -334,7 +336,7 @@ def PyExec(self):
         for ipk, peak in enumerate(peaks):
             prog_reporter.report("Integrating")
             peak_intens, peak_sigma = 0.0, 0.0
-            status = PEAK_STATUS.NO_PEAK
+            peak_status = PEAK_STATUS.NO_PEAK
 
             detid = peak.getDetectorID()
             bank_name = peaks.column("BankName")[ipk]
@@ -375,7 +377,7 @@ def PyExec(self):
                 continue  # skip peak
 
             # update peak mask based on I/sig from fit
-            *_, i_over_sigma = calc_intens_and_sigma_arrays(fit_result, error_strategy)
+            *_, i_over_sigma, _ = calc_intens_and_sigma_arrays(fit_result, error_strategy)
             non_bg_mask = np.zeros(peak_data.detids.shape, dtype=bool)
             non_bg_mask.flat[initial_peak_mask] = i_over_sigma > i_over_sig_threshold
             peak_mask = find_peak_cluster_in_window(non_bg_mask, (peak_data.irow, peak_data.icol))
@@ -384,7 +386,7 @@ def PyExec(self):
 
             is_on_edge = np.any(np.logical_and(peak_mask, peak_data.det_edges))
             if is_on_edge:
-                status = PEAK_STATUS.ON_EDGE
+                peak_status = PEAK_STATUS.ON_EDGE
                 if not integrate_on_edge:
                     self.delete_fit_result_workspaces(fit_result)
                     continue  # skip peak
@@ -406,19 +408,21 @@ def PyExec(self):
             peak_mask.flat[initial_peak_mask] = fit_mask
 
             # calculate intensity
-            status = PEAK_STATUS.VALID
-            intens, sigma, _ = calc_intens_and_sigma_arrays(fit_result, error_strategy)
+            peak_status = PEAK_STATUS.VALID
+            intens, sigma, _, peak_limits = calc_intens_and_sigma_arrays(fit_result, error_strategy)
             peak_intens = np.sum(intens[fit_mask])
             peak_sigma = np.sqrt(np.sum(sigma[fit_mask] ** 2))
 
+            self._set_peak_shape(peak, peak_data.detids[peak_mask], peak_limits[fit_mask])
+
             if output_file:
                 intens_over_sig = peak_intens / peak_sigma if peak_sigma > 0 else 0.0
                 results.append(
                     LineProfileResult(
                         ipk,
                         peak,
                         intens_over_sig,
-                        status,
+                        peak_status,
                         peak_mask,
                         fit_mask,
                         func_generator.ysum,
@@ -439,6 +443,15 @@ def PyExec(self):
         # assign output
         self.setProperty("OutputWorkspace", peaks)
 
+    def _set_peak_shape(self, peak, det_ids, fit_limits):
+        det_bin_list = []
+        for detid, limits in zip(det_ids, fit_limits):
+            if limits:
+                det_bin_list.append((int(detid), limits[0], limits[-1]))
+        if len(det_bin_list) > 0:
+            peak_shape = PeakShapeDetectorBin(det_bin_list, SpecialCoordinateSystem.NONE, self.name(), self.version())
+            peak.setPeakShape(peak_shape)
+
     def exec_child_alg(self, alg_name, **kwargs):
         alg = self.createChildAlgorithm(alg_name, enableLogging=False)
         alg.initialize()
@@ -478,7 +491,7 @@ def __init__(self, peak_params_to_fix: Sequence[str]):
         self.cen_par_name: str = None
         self.intens_par_name: str = None
         self.width_par_name: str = None
-        self.width_max: float = None
+        self.width_min: float = None
         self.width_max: float = None
         self.peak_params_to_fix: Sequence[str] = peak_params_to_fix
         self.peak_mask: np.ndarray[float] = None
@@ -732,7 +745,7 @@ def calc_sigma_from_summation(xdat, edat_sq, ypeak, cutoff=0.025):
     ilo = np.clip(np.argmin(abs(ypeak_cumsum - cutoff)), a_min=0, a_max=nbins // 2)
     ihi = np.clip(np.argmin(abs(ypeak_cumsum - (1 - cutoff))), a_min=nbins // 2, a_max=nbins - 1) + 1
     bin_width = np.diff(xdat[ilo : ihi + 1])
-    return np.sqrt(np.sum(edat_sq[ilo:ihi] * (bin_width**2)))
+    return np.sqrt(np.sum(edat_sq[ilo:ihi] * (bin_width**2))), (xdat[ilo], xdat[ihi])
 
 
 def calc_intens_and_sigma_arrays(fit_result, error_strategy):
@@ -741,6 +754,7 @@ def calc_intens_and_sigma_arrays(fit_result, error_strategy):
     sigma = np.zeros(intens.shape)
     intens_over_sig = np.zeros(intens.shape)
     peak_func = FunctionFactory.Instance().createPeakFunction(function[0][0].name())
+    peak_limits = np.full(intens.shape, None)
     for idom, comp_func in enumerate(function):
         [peak_func.setParameter(iparam, comp_func.getParameterValue(iparam)) for iparam in range(peak_func.nParams())]
         intens[idom] = peak_func.intensity()
@@ -749,10 +763,10 @@ def calc_intens_and_sigma_arrays(fit_result, error_strategy):
             sigma[idom] = peak_func.intensityError()
         else:
             ws_fit = get_eval_ws(fit_result["OutputWorkspace"], idom)
-            sigma[idom] = calc_sigma_from_summation(ws_fit.readX(0), ws_fit.readE(0) ** 2, ws_fit.readY(3))
+            sigma[idom], peak_limits[idom] = calc_sigma_from_summation(ws_fit.readX(0), ws_fit.readE(0) ** 2, ws_fit.readY(3))
     ivalid = ~np.isclose(sigma, 0)
     intens_over_sig[ivalid] = intens[ivalid] / sigma[ivalid]
-    return intens, sigma, intens_over_sig
+    return intens, sigma, intens_over_sig, peak_limits
 
 
 def get_eval_ws(out_ws_name, idom):

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

@@ -5,6 +5,9 @@
 #   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
 # SPDX - License - Identifier: GPL - 3.0 +
 import unittest
+from unittest import mock
+from unittest.mock import patch
+
 from mantid.simpleapi import (
     IntegratePeaks1DProfile,
     CreatePeaksWorkspace,
@@ -20,6 +23,7 @@
 import tempfile
 import shutil
 from os import path
+import json
 
 
 class IntegratePeaks1DProfileTest(unittest.TestCase):
@@ -86,6 +90,35 @@ def test_exec_IntegrateIfOnEdge_True(self):
         self.assertAlmostEqual(out.column("Intens/SigInt")[0], self.default_intens_over_sigma, delta=1e-1)
         self.assertAlmostEqual(out.column("Intens/SigInt")[1], self.default_intens_over_sigma, delta=1e-1)
 
+    def test_peak_shapes_when_fit_successful(self):
+        out = IntegratePeaks1DProfile(
+            InputWorkspace=self.ws, PeaksWorkspace=self.peaks_edge, OutputWorkspace="peaks_int_1", **self.profile_kwargs
+        )
+        self.assertEqual(len(out), 2)
+        for pk in out:
+            self.assertEqual(pk.getPeakShape().shapeName(), "detectorbin")
+            self.assertEqual(pk.getPeakShape().algorithmName(), "IntegratePeaks1DProfile")
+            self.assertEqual(pk.getPeakShape().algorithmVersion(), 1)
+            pk_shape_dict = json.loads(pk.getPeakShape().toJSON())
+            self.assertEqual(len(pk_shape_dict["detectors"]), 2)
+            for det in pk_shape_dict["detectors"]:
+                self.assertTrue(det["startX"] < det["endX"])
+
+    @patch("plugins.algorithms.IntegratePeaks1DProfile.IntegratePeaks1DProfile.delete_fit_result_workspaces")
+    @patch(
+        "plugins.algorithms.IntegratePeaks1DProfile.IntegratePeaks1DProfile.exec_fit",
+    )
+    def test_peak_shapes_when_fit_failed(self, mock_fit, mock_fit_result_workspaces):
+        mock_fit_return = {"success": False}
+        mock_fit.return_value = mock_fit_return
+        out = IntegratePeaks1DProfile(
+            InputWorkspace=self.ws, PeaksWorkspace=self.peaks_edge, OutputWorkspace="peaks_int_1", **self.profile_kwargs
+        )
+        self.assertEqual(len(out), 2)
+        for pk in out:
+            self.assertEqual(pk.getPeakShape().shapeName(), "none")
+        mock_fit_result_workspaces.assert_has_calls([mock.call(mock_fit_return), mock.call(mock_fit_return)])
+
     def test_exec_IntegrateIfOnEdge_False(self):
         kwargs = self.profile_kwargs.copy()
         kwargs["IntegrateIfOnEdge"] = False

docs/source/release/v6.13.0/Diffraction/Single_Crystal/New_features/38969.rst

@@ -0,0 +1,2 @@
+- Added ``detectorbin`` peak shape for the peaks integrated with :ref:`IntegratePeaks1DProfile <algm-IntegratePeaks1DProfile>` integration algorithm.
+- By accessing the detectorbin peak shape, users can now view the detector IDs and the corresponding range in the X dimension associated with each detector for each successfully integrated peak from the algorithm.