Merge pull request mantidproject#36481 from mantidproject/36322_shoeb…

…ox_integration_algorithm

New shoebox integration algorithm for single-crystal Bragg peaks

Framework/PythonInterface/plugins/CMakeLists.txt

@@ -93,6 +93,7 @@ set(PYTHON_PLUGINS
     algorithms/InelasticEMUauReduction.py
     algorithms/IntegratePeaksProfileFitting.py
     algorithms/IntegratePeaksSkew.py
+    algorithms/IntegratePeaksShoeboxTOF.py
     algorithms/LRAutoReduction.py
     algorithms/LRDirectBeamSort.py
     algorithms/LRPeakSelection.py

Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py

@@ -155,22 +155,14 @@ def PyExec(self):
             else:
                 nbins = self.getProperty("NBins").value
 
-            # make a kernel with background subtraction shell of approx. same number of elements as peak region
-            nrows_bg = max(1, nrows // 8)  # padding either side of e.g. nrows for bg shell
-            ncols_bg = max(1, ncols // 8)
-            nbins_bg = max(1, nbins // 8)
-            kernel = np.zeros((nrows + 2 * nrows_bg, ncols + 2 * ncols_bg, nbins + 2 * nbins_bg))
-            kernel[nrows_bg:-nrows_bg, ncols_bg:-ncols_bg, nbins_bg:-nbins_bg] = 1
-            bg_mask = np.logical_not(kernel)
-            kernel[bg_mask] = -(nrows * ncols * nbins) / bg_mask.sum()  # such that kernel.sum() = 0
-
             # get data in detector coords
             peak_data = array_converter.get_peak_data(dummy_pk, detid, bank.getName(), bank.xpixels(), bank.ypixels(), 1, 1)
-            _, y, e = peak_data.get_data_arrays()  # 3d arrays [rows x cols x tof]
+            _, y, esq, _ = peak_data.get_data_arrays()  # 3d arrays [rows x cols x tof]
             # perform convolutions to integrate kernel/shoebox
             # pad with nearest so don't get peaks at edge when -ve values go outside data extent
+            kernel = make_kernel(nrows, ncols, nbins)
             yconv = convolve(input=y, weights=kernel, mode="nearest")
-            econv = np.sqrt(convolve(input=e**2, weights=kernel**2, mode="nearest"))
+            econv = np.sqrt(convolve(input=esq, weights=kernel**2, mode="nearest"))
             with np.errstate(divide="ignore", invalid="ignore"):
                 intens_over_sig = yconv / econv  # ignore 0/0 which produces NaN (recall NaN > x = False)
 
@@ -180,8 +172,8 @@ def PyExec(self):
             labels, nlabels = label(pk_mask)  # identify contiguous nearest-neighbour connected regions
             # identify labels of peaks above min size
             min_size = int(min_frac_size * kernel.size)
-            nbins = sum_labels(pk_mask, labels, range(1, nlabels + 1))
-            ilabels = np.flatnonzero(nbins > min_size) + 1
+            npixels = sum_labels(pk_mask, labels, range(1, nlabels + 1))
+            ilabels = np.flatnonzero(npixels > min_size) + 1
 
             # find index of maximum in I/sigma for each valid peak (label index in ilabels)
             imaxs = maximum_position(intens_over_sig, labels, ilabels)
@@ -208,8 +200,11 @@ def PyExec(self):
                 # remove peaks on edge if required
                 do_add_peak = True
                 if remove_on_edge:
+                    nrows_bg = kernel.shape[0] - nrows
+                    ncols_bg = kernel.shape[1] - ncols
+                    nbins_bg = kernel.shape[2] - nbins
                     for idim, (index, nbg) in enumerate([(irow_max, nrows_bg), (icol_max, ncols_bg), (itof_max, nbins_bg)]):
-                        if not 2 * nbg <= index <= y.shape[idim] - 2 * nbg - 1:
+                        if not nbg <= index <= y.shape[idim] - nbg - 1:
                             do_add_peak = False
                             break
                 if do_add_peak:
@@ -247,5 +242,22 @@ def exec_child_alg(self, alg_name, **kwargs):
             return None
 
 
+def make_kernel(nrows, ncols, nbins):
+    # make a kernel with background subtraction shell of approx. same number of elements as peak region
+    kernel_shape, (nrows_bg, ncols_bg, nbins_bg) = get_kernel_shape(nrows, ncols, nbins)
+    kernel = np.zeros(kernel_shape)
+    kernel[nrows_bg:-nrows_bg, ncols_bg:-ncols_bg, nbins_bg:-nbins_bg] = 1
+    bg_mask = np.logical_not(kernel)
+    kernel[bg_mask] = -(nrows * ncols * nbins) / bg_mask.sum()  # such that kernel.sum() = 0
+    return kernel
+
+
+def get_kernel_shape(nrows, ncols, nbins):
+    nrows_bg = max(1, nrows // 16)  # padding either side of e.g. nrows for bg shell
+    ncols_bg = max(1, ncols // 16)
+    nbins_bg = max(1, nbins // 16)
+    return (nrows + 2 * nrows_bg, ncols + 2 * ncols_bg, nbins + 2 * nbins_bg), (nrows_bg, ncols_bg, nbins_bg)
+
+
 # register algorithm with mantid
 AlgorithmFactory.subscribe(FindSXPeaksConvolve)