From 9a90dd6e2ae5d878455b62ab333a1d6f6ed31091 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 6 Nov 2023 16:39:25 +0000 Subject: [PATCH 01/37] Move making kernel to separate function Also fixed bg shell calculation to make it closer to having same number elements as peak extent in kernel --- .../plugins/algorithms/FindSXPeaksConvolve.py | 27 ++++++++++++------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py index d2dd23e8698e..9eb9b6baa150 100644 --- a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py +++ b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py @@ -155,20 +155,12 @@ 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] # 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")) with np.errstate(divide="ignore", invalid="ignore"): @@ -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[0] - 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,17 @@ 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 + 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) + 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 + return kernel + + # register algorithm with mantid AlgorithmFactory.subscribe(FindSXPeaksConvolve) From 1542bbcd3b2cb2bdf5fece4ade8e476a580b2bee Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 8 Nov 2023 15:21:31 +0000 Subject: [PATCH 02/37] Return errors squared and ws index from PeakData.get_data_arrays --- .../PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py | 4 ++-- .../PythonInterface/plugins/algorithms/IntegratePeaksSkew.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py index 9eb9b6baa150..6006d9e45a06 100644 --- a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py +++ b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py @@ -157,12 +157,12 @@ def PyExec(self): # 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) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py index 80103770dbd6..e81109e9106a 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py @@ -266,7 +266,7 @@ def get_data_arrays(self): if len(xvals) > signal.shape[-1]: xvals = 0.5 * (xvals[:-1] + xvals[1:]) # convert to bin centers xcens[irow, icol, :] = xvals - return xcens, signal, errors + return xcens, signal, errors**2, ispecs def get_roi_on_detector(self, detids): isort = np.argsort(detids.flatten()) From c2a3b499b363e69bec316645d148c0bba57e6975 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 8 Nov 2023 15:28:22 +0000 Subject: [PATCH 03/37] Fix bug in FindSCPeaksConvolve introduced with make_kernel func --- .../plugins/algorithms/FindSXPeaksConvolve.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py index 6006d9e45a06..48ef08c74049 100644 --- a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py +++ b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py @@ -172,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) @@ -202,7 +202,7 @@ def PyExec(self): if remove_on_edge: nrows_bg = kernel.shape[0] - nrows ncols_bg = kernel.shape[1] - ncols - nbins_bg = kernel.shape[0] - nbins + 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 nbg <= index <= y.shape[idim] - nbg - 1: do_add_peak = False From 18c4a013fc5ddca6fc1c90c45a604b9dafc78182 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Nov 2023 12:10:21 +0000 Subject: [PATCH 04/37] Add get_kernel_shape function --- .../plugins/algorithms/FindSXPeaksConvolve.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py index 48ef08c74049..2fb3d1666677 100644 --- a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py +++ b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py @@ -244,15 +244,20 @@ def exec_child_alg(self, alg_name, **kwargs): def make_kernel(nrows, ncols, nbins): # make a kernel with background subtraction shell of approx. same number of elements as peak region - 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) - kernel = np.zeros((nrows + 2 * nrows_bg, ncols + 2 * ncols_bg, nbins + 2 * nbins_bg)) + 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) From 211c65b4b31a4897edf74e57065e65baf5af0bd4 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Nov 2023 16:37:04 +0000 Subject: [PATCH 05/37] Add new shoebox integration algorithm --- .../PythonInterface/plugins/CMakeLists.txt | 1 + .../algorithms/IntegratePeaksShoeboxTOF.py | 478 ++++++++++++++++++ 2 files changed, 479 insertions(+) create mode 100644 Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py diff --git a/Framework/PythonInterface/plugins/CMakeLists.txt b/Framework/PythonInterface/plugins/CMakeLists.txt index abd30b0ef3db..456f91455e7a 100644 --- a/Framework/PythonInterface/plugins/CMakeLists.txt +++ b/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 diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py new file mode 100644 index 000000000000..b2e5562d86d2 --- /dev/null +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -0,0 +1,478 @@ +# 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 +# SPDX - License - Identifier: GPL - 3.0 + +from mantid.api import ( + DataProcessorAlgorithm, + AlgorithmFactory, + MatrixWorkspaceProperty, + IPeaksWorkspaceProperty, + WorkspaceUnitValidator, + FileProperty, + FileAction, +) +from mantid.kernel import ( + Direction, + FloatBoundedValidator, + IntBoundedValidator, + StringListValidator, + EnabledWhenProperty, + PropertyCriterion, +) +import numpy as np +from scipy.ndimage import convolve +from IntegratePeaksSkew import InstrumentArrayConverter, get_fwhm_from_back_to_back_params +from FindSXPeaksConvolve import make_kernel, get_kernel_shape +from enum import Enum + + +class PEAK_STATUS(Enum): + WEAK = "Weak peak" + STRONG = "Strong peak" + ON_EDGE = "Peak mask is on the detector edge" + NO_PEAK = "No peak detected." + + +class ShoeboxResult: + """ + This class is used to hold data and integration parameters for single-crystal Bragg peaks + """ + + def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, status): + self.peak_shape = list(peak_shape) + self.kernel_shape = list(get_kernel_shape(*self.peak_shape)[0]) + self.ipos = list(ipos) + self.ipos_pk = list(ipk_pos) + self.labels = ["Row", "Col", "TOF"] + self.ysum = [] + self.extents = [] + # integrate y over each dim + for idim in range(len(peak_shape)): + self.ysum.append(y.sum(axis=idim)) + if idim < 2: + self.extents.append([0, y.shape[idim]]) + else: + self.extents.append([x[0], x[-1]]) + # extract peak properties + intens_over_sig = np.round(pk.getIntensityOverSigma(), 1) + hkl = np.round(pk.getHKL(), 2) + wl = np.round(pk.getWavelength(), 2) + tth = np.round(np.degrees(pk.getScattering()), 1) + d = np.round(pk.getDSpacing(), 2) + self.title = ( + f"{ipk} ({','.join(str(hkl)[1:-1].split())})" + f"\n$I/\\sigma$={intens_over_sig}\n" + rf"$\lambda$={wl} $\AA$; " + rf"$2\theta={tth}^\circ$; " + rf"d={d} $\AA$" + f"\n{status.value}" + ) + + def plot_integrated_peak(self, fig, axes, norm_func, rect_func): + for iax, ax in enumerate(axes): + extents = self.extents.copy() + extents.pop(iax) + im = ax.imshow( + self.ysum[iax], aspect=self.ysum[iax].shape[1] / self.ysum[iax].shape[0] + ) # , extent=np.flip(extents, axis=0).flatten()) + im.set_norm(norm_func()) + # plot shoebox + self.plot_shoebox(iax, ax, rect_func) + # plot peak position + cen = self.ipos_pk.copy() + cen.pop(iax) + ax.plot(cen[1], cen[0], "xr") + # set labels + labels = self.labels.copy() + labels.pop(iax) + ax.set_xlabel(labels[1]) + ax.set_ylabel(labels[0]) + # set title + fig.suptitle(self.title) + + def plot_shoebox(self, iax, ax, rect_func): + cen = self.ipos.copy() + cen.pop(iax) + # plot peak region + shape = self.peak_shape.copy() + shape.pop(iax) + self._plot_rect(ax, cen, shape, rect_func) + # plot kernel (incl. bg shell) + shape = self.kernel_shape.copy() + shape.pop(iax) + self._plot_rect(ax, cen, shape, rect_func, ls="--") + + def _plot_rect(self, ax, cen, shape, rect_func, ls="-"): + bottom_left = [cen[1] - shape[1] // 2 - 0.5, cen[0] - shape[0] // 2 - 0.5] + rect = rect_func(bottom_left, shape[1], shape[0], linewidth=1.5, edgecolor="r", facecolor="none", alpha=0.75, ls=ls) + ax.add_patch(rect) + + +class IntegratePeaksShoeboxTOF(DataProcessorAlgorithm): + def name(self): + return "IntegratePeaksShoeboxTOF" + + def category(self): + return "Diffraction\\Reduction" + + def seeAlso(self): + return ["IntegratePeaksSkew", "FindSXPeaksConvolve"] + + def summary(self): + return "Integrate single-crystal Bragg peaks in MatrixWorkspaces with x-unit of TOF using a shoebox." + + def PyInit(self): + # Input + self.declareProperty( + MatrixWorkspaceProperty( + name="InputWorkspace", defaultValue="", direction=Direction.Input, validator=WorkspaceUnitValidator("TOF") + ), + doc="A MatrixWorkspace to integrate (x-axis must be TOF).", + ) + self.declareProperty( + IPeaksWorkspaceProperty(name="PeaksWorkspace", defaultValue="", direction=Direction.Input), + doc="A PeaksWorkspace containing the peaks to integrate.", + ) + self.declareProperty( + IPeaksWorkspaceProperty(name="OutputWorkspace", defaultValue="", direction=Direction.Output), + doc="The output PeaksWorkspace will be a copy of the input PeaksWorkspace with the" " integrated intensities.", + ) + # shoebox dimensions + self.declareProperty( + name="NRows", + defaultValue=5, + direction=Direction.Input, + validator=IntBoundedValidator(lower=2), + doc="Number of row components in the detector to use in the convolution kernel. " + "For WISH row components correspond to pixels along a single tube.", + ) + self.declareProperty( + name="NCols", + defaultValue=5, + direction=Direction.Input, + validator=IntBoundedValidator(lower=2), + doc="Number of column components in the detector to use in the convolution kernel. " + "For WISH column components correspond to tubes.", + ) + self.declareProperty( + name="NBins", + defaultValue=10, + direction=Direction.Input, + validator=IntBoundedValidator(lower=2), + doc="Number of TOF bins to use in the convolution kernel.", + ) + self.declareProperty( + name="GetNBinsFromBackToBackParams", + defaultValue=False, + direction=Direction.Input, + doc="If true the number of TOF bins used in the convolution kernel will be calculated from the FWHM of the " + "BackToBackExponential peak using parameters defined in the instrument parameters.xml file.", + ) + self.declareProperty( + name="NFWHM", + defaultValue=4, + direction=Direction.Input, + validator=IntBoundedValidator(lower=1), + doc="If GetNBinsFromBackToBackParams=True then the number of TOF bins will be NFWHM x FWHM of the " + "BackToBackExponential at the peak detector and TOF.", + ) + use_nBins = EnabledWhenProperty("GetNBinsFromBackToBackParams", PropertyCriterion.IsDefault) + self.setPropertySettings("NBins", use_nBins) + use_nfwhm = EnabledWhenProperty("GetNBinsFromBackToBackParams", PropertyCriterion.IsNotDefault) + self.setPropertySettings("NFWHM", use_nfwhm) + self.declareProperty( + name="NShoeboxInWindow", + defaultValue=3.0, + direction=Direction.Input, + validator=FloatBoundedValidator(lower=2.0), + doc="Extent of window in TOF and detector row and column as a multiple of the shoebox length along each dimension.", + ) + # shoebox optimisation + self.declareProperty( + name="OptimiseShoebox", + defaultValue=True, + direction=Direction.Input, + doc="If OptimiseShoebox=True then shoebox size will be optimised to maximise Intensity/Sigma of the peak.", + ) + self.declareProperty( + name="WeakPeakStrategy", + defaultValue="Fix", + direction=Direction.Input, + validator=StringListValidator(["Fix", "NearestStrongPeak"]), + doc="If WeakPeakStrategy=Fix then fix the shoebox dimensions. If WeakPeakStrategy=NearestStrongPeak then" + "the shoebox dimensions will be taken from the nearest strong peak (determined by StrongPeakThreshold)." + "The TOF extent of the shoebox will be scaled by the FWHM of the BackToBackExponential peaks if" + "GetNBinsFromBackToBackParams=True.", + ) + self.declareProperty( + name="WeakPeakThreshold", + defaultValue=10.0, + direction=Direction.Input, + validator=FloatBoundedValidator(lower=0.0), + doc="Intenisty/Sigma threshold below which a peak is considered weak.", + ) + # peak validation + self.declareProperty( + name="IntegrateIfOnEdge", + defaultValue=False, + direction=Direction.Input, + doc="If IntegrateIfOnEdge=False then peaks with shoebox that includes detector IDs at the edge of the bank " + " will not be integrated.", + ) + self.declareProperty( + name="NRowsEdge", + defaultValue=1, + direction=Direction.Input, + validator=IntBoundedValidator(lower=1), + doc="Shoeboxes containing detectors NRowsEdge from the detector edge are defined as on the edge.", + ) + self.declareProperty( + name="NColsEdge", + defaultValue=1, + direction=Direction.Input, + validator=IntBoundedValidator(lower=1), + doc="Shoeboxes containing detectors NColsEdge from the detector edge are defined as on the edge.", + ) + # plotting + self.declareProperty( + FileProperty("OutputFile", "", FileAction.OptionalSave, ".pdf"), + "Optional file path in which to write diagnostic plots (note this will slow the " "execution of algorithm).", + ) + # Corrections + self.declareProperty( + name="LorentzCorrection", + defaultValue=True, + direction=Direction.Input, + doc="Correct the integrated intensity by multiplying by the Lorentz factor " + "sin(theta)^2 / lambda^4 - do not do this if the data have already been corrected.", + ) + self.setPropertyGroup("LorentzCorrection", "Corrections") + + def PyExec(self): + # get input + ws = self.getProperty("InputWorkspace").value + peaks = self.getProperty("PeaksWorkspace").value + # shoebox dimensions + nrows = self.getProperty("NRows").value + ncols = self.getProperty("NCols").value + get_nbins_from_b2bexp_params = self.getProperty("GetNBinsFromBackToBackParams").value + nfwhm = self.getProperty("NFWHM").value + nshoebox = self.getProperty("NShoeboxInWindow").value + # shoebox optimisation + do_optimise_shoebox = self.getProperty("OptimiseShoebox").value + weak_peak_strategy = self.getProperty("WeakPeakStrategy").value + weak_peak_threshold = self.getProperty("WeakPeakThreshold").value + # validation + integrate_on_edge = self.getProperty("IntegrateIfOnEdge").value + nrows_edge = self.getProperty("NRowsEdge").value + ncols_edge = self.getProperty("NColsEdge").value + # corrections + do_lorz_cor = self.getProperty("LorentzCorrection").value + # saving file + output_file = self.getProperty("OutputFile").value + + # create output table workspace + peaks = self.exec_child_alg("CloneWorkspace", InputWorkspace=peaks, OutputWorkspace="out_peaks") + + array_converter = InstrumentArrayConverter(ws) + ipk_weak = [] + results = np.empty(peaks.getNumberPeaks(), dtype=ShoeboxResult) + for ipk, peak in enumerate(peaks): + detid = peak.getDetectorID() + bank_name = peaks.column("BankName")[ipk] + pk_tof = peak.getTOF() + + # get shoebox kernel for initial integration + ispec = ws.getIndicesFromDetectorIDs([detid])[0] + itof = ws.binIndexOf(pk_tof, ispec) + bin_width = np.diff(ws.readX(ispec)[itof : itof + 2])[0] # used later to scale intensity + if get_nbins_from_b2bexp_params: + fwhm = get_fwhm_from_back_to_back_params(peak, ws, detid) + + nbins = max(3, int(nfwhm * fwhm / bin_width)) if fwhm is not None else self.getProperty("NBins").value + if not nbins % 2: + nbins += 1 # force to be odd + else: + nbins = self.getProperty("NBins").value + kernel = make_kernel(nrows, ncols, nbins) + + # get data array and crop + peak_data = array_converter.get_peak_data( + peak, detid, bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge + ) + x, y, esq, ispecs = peak_data.get_data_arrays() # 3d arrays [rows x cols x tof] + x = x[peak_data.irow, peak_data.icol, :] # take x at peak centre, should be same for all detectors + + # crop data array to TOF region of peak using shoebox dimension + itof = ws.yIndexOfX(pk_tof, ispec) # need index in y now (note x values are points even if were edges) + tof_slice = slice( + int(np.clip(itof - nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), + int(np.clip(itof + nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), + ) + + x = x[tof_slice] + y = y[:, :, tof_slice] + esq = esq[:, :, tof_slice] + + # perform initial integration + intens_over_sig = convolve_shoebox(y, esq, kernel) + + # identify best shoebox position near peak + ix = np.argmin(abs(x - pk_tof)) + ipos = find_nearest_peak_in_data(intens_over_sig, ispecs, x, ws, peaks, ipk, peak_data.irow, peak_data.icol, ix) + + # perform final integration if required + det_edges = peak_data.det_edges if not integrate_on_edge else None + if ipos is None: + status = PEAK_STATUS.NO_PEAK + intens, sigma = 0.0, 0.0 + else: + # integrate at that position (without smoothing I/sigma) + intens, sigma = integrate_shoebox_at_pos(y, esq, kernel, ipos, det_edges) + if np.isfinite(sigma): + status = PEAK_STATUS.STRONG if intens / sigma > weak_peak_threshold else PEAK_STATUS.WEAK + if status == PEAK_STATUS.STRONG and do_optimise_shoebox: + kernel, (nrows, ncols, nbins) = optimise_shoebox(y, esq, kernel.shape, ipos) + # re-integrate but this time check for overlap with edge + intens, sigma = integrate_shoebox_at_pos(y, esq, kernel, ipos, det_edges) + else: + status = PEAK_STATUS.ON_EDGE + intens, sigma = 0.0, 0.0 + + if status == PEAK_STATUS.WEAK and weak_peak_strategy == "NearestStrongPeak": + ipk_weak.append(ipk) # store pk index (and pos FWHM?) + else: + set_peak_intensity(peak, intens, sigma, bin_width, do_lorz_cor) + if output_file: + # save result for plotting + results[ipk] = ShoeboxResult(ipk, peak, x, y, [nrows, ncols, nbins], ipos, [peak_data.irow, peak_data.icol, ix], status) + # loop over weak peaks + # look for peaks with detectors in ispec window + # if no peaks found + # look for nearest strong peak + + # plot output + if output_file: + from matplotlib.pyplot import subplots, close + from matplotlib.patches import Rectangle + from matplotlib.colors import LogNorm + from matplotlib.backends.backend_pdf import PdfPages + + try: + with PdfPages(output_file) as pdf: + for result in results: + fig, axes = subplots(1, 3, figsize=(12, 5), subplot_kw={"projection": "mantid"}) + fig.subplots_adjust(wspace=0.3) # ensure plenty space between subplots (want to avoid slow tight_layout) + result.plot_integrated_peak(fig, axes, LogNorm, Rectangle) + pdf.savefig(fig) + close(fig) + except OSError: + raise RuntimeError( + f"OutputFile ({output_file}) could not be opened - please check it is not open by " + f"another programme and that the user has permission to write to that directory." + ) + + # assign output + self.setProperty("OutputWorkspace", peaks) + + def exec_child_alg(self, alg_name, **kwargs): + alg = self.createChildAlgorithm(alg_name, enableLogging=False) + alg.initialize() + for prop, value in kwargs.items(): + alg.setProperty(prop, value) + alg.execute() + if "OutputWorkspace" in alg.outputProperties(): + return alg.getProperty("OutputWorkspace").value + else: + return None + + +def convolve_shoebox(y, esq, kernel): + yconv = convolve(input=y, weights=kernel, 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 + intens_over_sig[~np.isfinite(intens_over_sig)] = 0 + intens_over_sig = convolve(intens_over_sig, weights=np.ones(3 * [3]) / (3**3), mode="constant") # 0 pad + return intens_over_sig + + +def find_nearest_peak_in_data(data, ispecs, x, ws, peaks, ipk, irow, icol, ix): + tofs = peaks.column("TOF") + ispecs_peaks = ws.getIndicesFromDetectorIDs([int(p.getDetectorID()) for p in peaks]) + ipks_near = np.flatnonzero(np.logical_and.reduce((tofs >= x.min(), tofs <= x.max(), np.isin(ispecs_peaks, ispecs)))) + ipks_near = np.delete(ipks_near, np.where(ipks_near == ipk)) # remove the peak of interest + if len(ipks_near) > 0: + # get position of nearby peaks in data array in fractional coords + shape = np.array(data.shape) + pos_near = [np.array(*np.where(ispecs == ispecs_peaks[ii]), np.array(np.argmin(x - tofs[ii]))) / shape for ii in ipks_near] + # sort data in descending order and select strongest peak nearest to pk (in fractional coordinates) + isort = list(zip(*np.unravel_index(np.argsort(-data, axis=None), data.shape))) + pk_pos = np.array([irow, icol, ix]) / np.array(data.shape) + imax_nearest = None + for ibin in isort: + # calc distance to pk pos + bin_pos = np.array(ibin) / shape + pk_dist_sq = np.sum(pk_pos - bin_pos) + for pos in pos_near: + if np.sum(pos - bin_pos) > pk_dist_sq: + imax_nearest = ibin + break + else: + continue # executed if inner loop did not break (i.e. pixel closer to nearby peak than this peak) + break # execute if inner loop did break and else branch ignored (i.e. found bin closest to this peak) + return imax_nearest # could be None if no peak found + else: + # no nearby peaks - return position of maximum in data + return np.unravel_index(np.argmax(data), data.shape) + + +def integrate_shoebox_at_pos(y, esq, kernel, ipos, det_edges=None): + slices = tuple([slice(ii - shape // 2, ii + shape // 2 + 1) for ii, shape in zip(ipos, kernel.shape)]) + intens, sigma = 2 * [np.nan] + if det_edges is None or not det_edges[slices[:-1]].any(): + intens = np.sum(y[slices] * kernel) + sigma = np.sqrt(np.sum(esq[slices] * (kernel**2))) + return intens, sigma + + +def optimise_shoebox(y, esq, current_shape, ipos): + # find limits of kernel size from data size + lengths = [] + for idim, length in enumerate(y.shape): + max_length = 2 * min(ipos[idim], length - ipos[idim]) + 1 + min_length = max(3, current_shape[idim] // 2) + if not min_length % 2: + min_length += 1 # force odd + lengths.append(np.arange(min_length, max_length - 2 * max(1, max_length // 16), 2)) + # loop over all possible kernel sizes + best_peak_shape = None + best_kernel = None + best_i_over_sig = -np.inf # null value + for nrows in lengths[0]: + for ncols in lengths[1]: + for nbins in lengths[2]: + kernel = make_kernel(nrows, ncols, nbins) + intens, sigma = integrate_shoebox_at_pos(y, esq, kernel, ipos) + i_over_sig = intens / sigma + if i_over_sig > best_i_over_sig: + best_i_over_sig = i_over_sig + best_peak_shape = (nrows, ncols, nbins) + best_kernel = kernel + return best_kernel, best_peak_shape + + +def set_peak_intensity(pk, intens, sigma, bin_width, do_lorz_cor): + if do_lorz_cor: + L = (np.sin(pk.getScattering() / 2) ** 2) / (pk.getWavelength() ** 4) # at updated peak pos + else: + L = 1 + # set peak object intensity + pk.setIntensity(L * intens * bin_width) + pk.setSigmaIntensity(L * sigma * bin_width) + + +# register algorithm with mantid +AlgorithmFactory.subscribe(IntegratePeaksShoeboxTOF) From 978f1b1f902ac100dd7a93c019abf5c263fd313b Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Nov 2023 17:04:26 +0000 Subject: [PATCH 06/37] Add additional input enabling and grouping in init --- .../algorithms/IntegratePeaksShoeboxTOF.py | 33 +++++++++++++------ 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index b2e5562d86d2..af2fd0d1276d 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -144,7 +144,7 @@ def PyInit(self): name="NRows", defaultValue=5, direction=Direction.Input, - validator=IntBoundedValidator(lower=2), + validator=IntBoundedValidator(lower=3), doc="Number of row components in the detector to use in the convolution kernel. " "For WISH row components correspond to pixels along a single tube.", ) @@ -152,15 +152,15 @@ def PyInit(self): name="NCols", defaultValue=5, direction=Direction.Input, - validator=IntBoundedValidator(lower=2), + validator=IntBoundedValidator(lower=3), doc="Number of column components in the detector to use in the convolution kernel. " "For WISH column components correspond to tubes.", ) self.declareProperty( name="NBins", - defaultValue=10, + defaultValue=11, direction=Direction.Input, - validator=IntBoundedValidator(lower=2), + validator=IntBoundedValidator(lower=3), doc="Number of TOF bins to use in the convolution kernel.", ) self.declareProperty( @@ -189,6 +189,8 @@ def PyInit(self): validator=FloatBoundedValidator(lower=2.0), doc="Extent of window in TOF and detector row and column as a multiple of the shoebox length along each dimension.", ) + for prop in ["NRows", "NCols", "NBins", "GetNBinsFromBackToBackParams", "NFWHM", "NShoeboxInWindow"]: + self.setPropertyGroup(prop, "Shoebox Dimensions") # shoebox optimisation self.declareProperty( name="OptimiseShoebox", @@ -208,11 +210,16 @@ def PyInit(self): ) self.declareProperty( name="WeakPeakThreshold", - defaultValue=10.0, + defaultValue=0.0, direction=Direction.Input, validator=FloatBoundedValidator(lower=0.0), doc="Intenisty/Sigma threshold below which a peak is considered weak.", ) + optimise_shoebox_enabled = EnabledWhenProperty("OptimiseShoebox", PropertyCriterion.IsDefault) + self.setPropertySettings("WeakPeakStrategy", optimise_shoebox_enabled) + self.setPropertySettings("WeakPeakThreshold", optimise_shoebox_enabled) + for prop in ["OptimiseShoebox", "WeakPeakStrategy", "WeakPeakThreshold"]: + self.setPropertyGroup(prop, "Shoebox Optimisation") # peak validation self.declareProperty( name="IntegrateIfOnEdge", @@ -235,11 +242,11 @@ def PyInit(self): validator=IntBoundedValidator(lower=1), doc="Shoeboxes containing detectors NColsEdge from the detector edge are defined as on the edge.", ) - # plotting - self.declareProperty( - FileProperty("OutputFile", "", FileAction.OptionalSave, ".pdf"), - "Optional file path in which to write diagnostic plots (note this will slow the " "execution of algorithm).", - ) + edge_check_enabled = EnabledWhenProperty("IntegrateIfOnEdge", PropertyCriterion.IsDefault) + self.setPropertySettings("NRowsEdge", edge_check_enabled) + self.setPropertySettings("NColsEdge", edge_check_enabled) + for prop in ["IntegrateIfOnEdge", "NRowsEdge", "NColsEdge"]: + self.setPropertyGroup(prop, "Edge Checking") # Corrections self.declareProperty( name="LorentzCorrection", @@ -249,6 +256,12 @@ def PyInit(self): "sin(theta)^2 / lambda^4 - do not do this if the data have already been corrected.", ) self.setPropertyGroup("LorentzCorrection", "Corrections") + # plotting + self.declareProperty( + FileProperty("OutputFile", "", FileAction.OptionalSave, ".pdf"), + "Optional file path in which to write diagnostic plots (note this will slow the " "execution of algorithm).", + ) + self.setPropertyGroup("OutputFile", "Plotting") def PyExec(self): # get input From 5753df3def7b5cf474d7df3ecac4ed37b4ffd3b1 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Nov 2023 17:07:08 +0000 Subject: [PATCH 07/37] Fix bug appending weak peak shoebox to be optimised --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index af2fd0d1276d..4a67346f291d 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -354,7 +354,7 @@ def PyExec(self): status = PEAK_STATUS.ON_EDGE intens, sigma = 0.0, 0.0 - if status == PEAK_STATUS.WEAK and weak_peak_strategy == "NearestStrongPeak": + if status == PEAK_STATUS.WEAK and do_optimise_shoebox and weak_peak_strategy == "NearestStrongPeak": ipk_weak.append(ipk) # store pk index (and pos FWHM?) else: set_peak_intensity(peak, intens, sigma, bin_width, do_lorz_cor) From 247e57dd4936fb13f99301818e6b792cfb780f51 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Nov 2023 17:16:05 +0000 Subject: [PATCH 08/37] Add validateInputs method for custom input validation --- .../algorithms/IntegratePeaksShoeboxTOF.py | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 4a67346f291d..30d00a799eee 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -263,6 +263,29 @@ def PyInit(self): ) self.setPropertyGroup("OutputFile", "Plotting") + def validateInputs(self): + issues = dict() + # check shoebox dimensions + for prop in ["NRows", "NCols", "NBins"]: + if not self.getProperty(prop).value % 2: + issues[prop] = f"{prop} must be an odd number." + # check valid peak workspace + ws = self.getProperty("InputWorkspace").value + inst = ws.getInstrument() + pk_ws = self.getProperty("PeaksWorkspace").value + if inst.getName() != pk_ws.getInstrument().getName(): + issues["PeaksWorkspace"] = "PeaksWorkspace must have same instrument as the InputWorkspace." + if pk_ws.getNumberPeaks() < 1: + issues["PeaksWorkspace"] = "PeaksWorkspace must have at least 1 peak." + # check that is getting dTOF from back-to-back params then they are present in instrument + if self.getProperty("GetNBinsFromBackToBackParams").value: + # check at least first peak in workspace has back to back params + if not inst.getComponentByName(pk_ws.column("BankName")[0]).hasParameter("B"): + issues[ + "GetNBinsFromBackToBackParams" + ] = "Workspace doesn't have back to back exponential coefficients defined in the parameters.xml file." + return issues + def PyExec(self): # get input ws = self.getProperty("InputWorkspace").value From 08e411b7e3459404bd9a38a5982077f89673f6c4 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Nov 2023 17:31:52 +0000 Subject: [PATCH 09/37] Fix bug with shoebox size overwritten if optimised for previous peak --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 30d00a799eee..11c7c871731e 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -291,8 +291,6 @@ def PyExec(self): ws = self.getProperty("InputWorkspace").value peaks = self.getProperty("PeaksWorkspace").value # shoebox dimensions - nrows = self.getProperty("NRows").value - ncols = self.getProperty("NCols").value get_nbins_from_b2bexp_params = self.getProperty("GetNBinsFromBackToBackParams").value nfwhm = self.getProperty("NFWHM").value nshoebox = self.getProperty("NShoeboxInWindow").value @@ -321,6 +319,8 @@ def PyExec(self): pk_tof = peak.getTOF() # get shoebox kernel for initial integration + nrows = self.getProperty("NRows").value # get these inside loop as overwritten if shoebox optimised + ncols = self.getProperty("NCols").value ispec = ws.getIndicesFromDetectorIDs([detid])[0] itof = ws.binIndexOf(pk_tof, ispec) bin_width = np.diff(ws.readX(ispec)[itof : itof + 2])[0] # used later to scale intensity @@ -347,7 +347,6 @@ def PyExec(self): int(np.clip(itof - nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), int(np.clip(itof + nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), ) - x = x[tof_slice] y = y[:, :, tof_slice] esq = esq[:, :, tof_slice] From 5ae13fba57d507b639a16467832e82b7880e2c8a Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 14 Nov 2023 10:15:48 +0000 Subject: [PATCH 10/37] Add functionality to find shoebox dims from nearest strong peak --- .../algorithms/IntegratePeaksShoeboxTOF.py | 232 +++++++++++++----- 1 file changed, 171 insertions(+), 61 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 11c7c871731e..256ed43b9c7f 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -35,12 +35,21 @@ class PEAK_STATUS(Enum): NO_PEAK = "No peak detected." +class WeakPeak: + def __init__(self, ipk, ispec, tof, tof_fwhm, ipks_near): + self.ipk = ipk + self.ispec = ispec # spectrum at center of kernel that corresponds to max I/sigma from convolution + self.tof = tof # TOF at center of kernel that corresponds to max I/sigma from convolution + self.tof_fwhm = tof_fwhm + self.ipks_near = ipks_near + + class ShoeboxResult: """ This class is used to hold data and integration parameters for single-crystal Bragg peaks """ - def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, status): + def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, status, strong_peak=None, ipk_strong=None): self.peak_shape = list(peak_shape) self.kernel_shape = list(get_kernel_shape(*self.peak_shape)[0]) self.ipos = list(ipos) @@ -48,6 +57,7 @@ def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, status): self.labels = ["Row", "Col", "TOF"] self.ysum = [] self.extents = [] + self.status = status # integrate y over each dim for idim in range(len(peak_shape)): self.ysum.append(y.sum(axis=idim)) @@ -69,6 +79,8 @@ def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, status): rf"d={d} $\AA$" f"\n{status.value}" ) + if strong_peak is not None and ipk_strong is not None: + self.title = f"{self.title} (shoebox taken from ipk={ipk_strong} {np.round(strong_peak.getHKL(), 2)})" def plot_integrated_peak(self, fig, axes, norm_func, rect_func): for iax, ax in enumerate(axes): @@ -311,8 +323,9 @@ def PyExec(self): peaks = self.exec_child_alg("CloneWorkspace", InputWorkspace=peaks, OutputWorkspace="out_peaks") array_converter = InstrumentArrayConverter(ws) - ipk_weak = [] - results = np.empty(peaks.getNumberPeaks(), dtype=ShoeboxResult) + weak_peaks_list = [] + ipks_strong = [] + results = np.full(peaks.getNumberPeaks(), None) for ipk, peak in enumerate(peaks): detid = peak.getDetectorID() bank_name = peaks.column("BankName")[ipk] @@ -328,8 +341,7 @@ def PyExec(self): fwhm = get_fwhm_from_back_to_back_params(peak, ws, detid) nbins = max(3, int(nfwhm * fwhm / bin_width)) if fwhm is not None else self.getProperty("NBins").value - if not nbins % 2: - nbins += 1 # force to be odd + nbins = round_up_to_odd_number(nbins) else: nbins = self.getProperty("NBins").value kernel = make_kernel(nrows, ncols, nbins) @@ -338,25 +350,14 @@ def PyExec(self): peak_data = array_converter.get_peak_data( peak, detid, bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge ) - x, y, esq, ispecs = peak_data.get_data_arrays() # 3d arrays [rows x cols x tof] - x = x[peak_data.irow, peak_data.icol, :] # take x at peak centre, should be same for all detectors - - # crop data array to TOF region of peak using shoebox dimension - itof = ws.yIndexOfX(pk_tof, ispec) # need index in y now (note x values are points even if were edges) - tof_slice = slice( - int(np.clip(itof - nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), - int(np.clip(itof + nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), - ) - x = x[tof_slice] - y = y[:, :, tof_slice] - esq = esq[:, :, tof_slice] + x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox) # perform initial integration intens_over_sig = convolve_shoebox(y, esq, kernel) # identify best shoebox position near peak ix = np.argmin(abs(x - pk_tof)) - ipos = find_nearest_peak_in_data(intens_over_sig, ispecs, x, ws, peaks, ipk, peak_data.irow, peak_data.icol, ix) + ipos = find_nearest_peak_in_data_window(intens_over_sig, ispecs, x, ws, peaks, ipk, peak_data.irow, peak_data.icol, ix) # perform final integration if required det_edges = peak_data.det_edges if not integrate_on_edge else None @@ -365,49 +366,94 @@ def PyExec(self): intens, sigma = 0.0, 0.0 else: # integrate at that position (without smoothing I/sigma) - intens, sigma = integrate_shoebox_at_pos(y, esq, kernel, ipos, det_edges) - if np.isfinite(sigma): - status = PEAK_STATUS.STRONG if intens / sigma > weak_peak_threshold else PEAK_STATUS.WEAK - if status == PEAK_STATUS.STRONG and do_optimise_shoebox: - kernel, (nrows, ncols, nbins) = optimise_shoebox(y, esq, kernel.shape, ipos) - # re-integrate but this time check for overlap with edge - intens, sigma = integrate_shoebox_at_pos(y, esq, kernel, ipos, det_edges) - else: - status = PEAK_STATUS.ON_EDGE - intens, sigma = 0.0, 0.0 + intens, sigma, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges) + if status == PEAK_STATUS.STRONG and do_optimise_shoebox: + kernel, (nrows, ncols, nbins) = optimise_shoebox(y, esq, kernel.shape, ipos) + # 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) if status == PEAK_STATUS.WEAK and do_optimise_shoebox and weak_peak_strategy == "NearestStrongPeak": - ipk_weak.append(ipk) # store pk index (and pos FWHM?) + # 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) + weak_peaks_list.append(WeakPeak(ipk, ispecs[ipos[0], ipos[1]], x[ipos[-1]], fwhm, ipks_near)) else: + if status == PEAK_STATUS.STRONG: + ipks_strong.append(ipk) set_peak_intensity(peak, intens, sigma, bin_width, do_lorz_cor) if output_file: # save result for plotting results[ipk] = ShoeboxResult(ipk, peak, x, y, [nrows, ncols, nbins], ipos, [peak_data.irow, peak_data.icol, ix], status) - # loop over weak peaks - # look for peaks with detectors in ispec window - # if no peaks found - # look for nearest strong peak - # plot output - if output_file: - from matplotlib.pyplot import subplots, close - from matplotlib.patches import Rectangle - from matplotlib.colors import LogNorm - from matplotlib.backends.backend_pdf import PdfPages - - try: - with PdfPages(output_file) as pdf: - for result in results: - fig, axes = subplots(1, 3, figsize=(12, 5), subplot_kw={"projection": "mantid"}) - fig.subplots_adjust(wspace=0.3) # ensure plenty space between subplots (want to avoid slow tight_layout) - result.plot_integrated_peak(fig, axes, LogNorm, Rectangle) - pdf.savefig(fig) - close(fig) - except OSError: - raise RuntimeError( - f"OutputFile ({output_file}) could not be opened - please check it is not open by " - f"another programme and that the user has permission to write to that directory." + if len(ipks_strong): + # set function for calculating distance metric between peaks + # if know back to back params then can scale TOF extent by ratio of FWHM + # otherwise just look for peak closest in QLab + calc_dist_func = calc_angle_between_peaks if get_nbins_from_b2bexp_params else calc_dQsq_between_peaks + + for weak_pk in weak_peaks_list: + # get peak + ipk = weak_pk.ipk + peak = peaks.getPeak(ipk) + bank_name = peaks.column("BankName")[ipk] + pk_tof = peak.getTOF() + # find nearest strong peak to get shoebox dimensions from + ipks_near_strong = [] + for ipk_near in weak_pk.ipks_near: + if results[ipk_near].status == PEAK_STATUS.STRONG: + ipks_near_strong.append(ipk_near) + if not ipks_near_strong: + # no peaks in detector window at any TOF, look in all table + ipks_near_strong = ipks_strong + # loop over strong peaks and find nearest peak + ipk_strong = None + dist_min = np.inf + for ipk_near in ipks_near_strong: + dist = calc_dist_func(peak, peaks.getPeak(ipk_near)) + if dist < dist_min: + ipk_strong = ipk_near + strong_pk = peaks.getPeak(ipk_strong) + + # get peak shape and make kernel + nrows, ncols, nbins = results[ipk_strong].peak_shape + if get_nbins_from_b2bexp_params: + # scale TOF extent by ratio of fwhm + strong_pk_fwhm = get_fwhm_from_back_to_back_params(strong_pk, ws, strong_pk.getDetectorID()) + nbins = max(3, int(nbins * (weak_pk.tof_fwhm / strong_pk_fwhm))) + nbins = round_up_to_odd_number(nbins) + kernel = make_kernel(nrows, ncols, nbins) + # get data array in peak region (keep same window size, nshoebox, for plotting) + peak_data = array_converter.get_peak_data( + peak, peak.getDetectorID(), bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge ) + x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox) + # 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) + set_peak_intensity(peak, intens, sigma, bin_width, do_lorz_cor) + if output_file: + # save result for plotting + results[ipk] = ShoeboxResult( + ipk, + peak, + x, + y, + [nrows, ncols, nbins], + ipos, + [peak_data.irow, peak_data.icol, np.argmin(abs(x - pk_tof))], + status, + strong_peak=strong_pk, + ipk_strong=ipk_strong, + ) + elif weak_peak_strategy == "NearestStrongPeak": + raise ValueError( + f"No peaks found with I/sigma > WeakPeakThreshold ({weak_peak_threshold}) - can't " + f"estimate shoebox dimension for weak peaks. Try reducing WeakPeakThreshold or set " + f"WeakPeakStrategy to Fix" + ) + + # plot output + plot_integration_reuslts(output_file, results) # assign output self.setProperty("OutputWorkspace", peaks) @@ -424,6 +470,58 @@ def exec_child_alg(self, alg_name, **kwargs): return None +def round_up_to_odd_number(number): + if not number % 2: + number += 1 + return number + + +def plot_integration_reuslts(output_file, results): + # import inside this function as not allowed to import at point algorithms are registered + from matplotlib.pyplot import subplots, close + from matplotlib.patches import Rectangle + from matplotlib.colors import LogNorm + from matplotlib.backends.backend_pdf import PdfPages + + try: + with PdfPages(output_file) as pdf: + for result in results: + if result: + fig, axes = subplots(1, 3, figsize=(12, 5), subplot_kw={"projection": "mantid"}) + fig.subplots_adjust(wspace=0.3) # ensure plenty space between subplots (want to avoid slow tight_layout) + result.plot_integrated_peak(fig, axes, LogNorm, Rectangle) + pdf.savefig(fig) + close(fig) + except OSError: + raise RuntimeError( + f"OutputFile ({output_file}) could not be opened - please check it is not open by " + f"another programme and that the user has permission to write to that directory." + ) + + +def calc_angle_between_peaks(pk1, pk2): + return abs(pk1.getQLabFrame().angle(pk2.getQLabFrame())) + + +def calc_dQsq_between_peaks(pk1, pk2): + return (pk1.getQLabFrame() - pk2.getQLabFrame()).norm2() + + +def get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox): + x, y, esq, ispecs = peak_data.get_data_arrays() # 3d arrays [rows x cols x tof] + x = x[peak_data.irow, peak_data.icol, :] # take x at peak centre, should be same for all detectors + # crop data array to TOF region of peak using shoebox dimension + itof = ws.yIndexOfX(pk_tof, ispec) # need index in y now (note x values are points even if were edges) + tof_slice = slice( + int(np.clip(itof - nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), + int(np.clip(itof + nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), + ) + x = x[tof_slice] + y = y[:, :, tof_slice] + esq = esq[:, :, tof_slice] + return x, y, esq, ispecs + + def convolve_shoebox(y, esq, kernel): yconv = convolve(input=y, weights=kernel, mode="nearest") econv = np.sqrt(convolve(input=esq, weights=kernel**2, mode="nearest")) @@ -434,12 +532,10 @@ def convolve_shoebox(y, esq, kernel): return intens_over_sig -def find_nearest_peak_in_data(data, ispecs, x, ws, peaks, ipk, irow, icol, ix): - tofs = peaks.column("TOF") - ispecs_peaks = ws.getIndicesFromDetectorIDs([int(p.getDetectorID()) for p in peaks]) - ipks_near = np.flatnonzero(np.logical_and.reduce((tofs >= x.min(), tofs <= x.max(), np.isin(ispecs_peaks, ispecs)))) - ipks_near = np.delete(ipks_near, np.where(ipks_near == ipk)) # remove the peak of interest +def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol, ix): + ipks_near, ispecs_peaks = find_ipks_in_window(ws, peaks, ispecs, ipk, tof_min=x.min(), tof_max=x.max()) if len(ipks_near) > 0: + tofs = peaks.column("TOF") # get position of nearby peaks in data array in fractional coords shape = np.array(data.shape) pos_near = [np.array(*np.where(ispecs == ispecs_peaks[ii]), np.array(np.argmin(x - tofs[ii]))) / shape for ii in ipks_near] @@ -464,13 +560,27 @@ def find_nearest_peak_in_data(data, ispecs, x, ws, peaks, ipk, irow, icol, ix): return np.unravel_index(np.argmax(data), data.shape) -def integrate_shoebox_at_pos(y, esq, kernel, ipos, det_edges=None): +def find_ipks_in_window(ws, peaks, ispecs, ipk, tof_min=None, tof_max=None): + ispecs_peaks = ws.getIndicesFromDetectorIDs([int(p.getDetectorID()) for p in peaks]) + ipks_near = np.isin(ispecs_peaks, ispecs) + if tof_min and tof_max: + tofs = peaks.column("TOF") + ipks_near = np.logical_and.reduce((tofs >= tof_min, tofs <= tof_max, ipks_near)) + ipks_near = np.flatnonzero(ipks_near) # convert to index + ipks_near = np.delete(ipks_near, np.where(ipks_near == ipk)) # remove the peak of interest + return ipks_near, ispecs_peaks + + +def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges=None): slices = tuple([slice(ii - shape // 2, ii + shape // 2 + 1) for ii, shape in zip(ipos, kernel.shape)]) - intens, sigma = 2 * [np.nan] if det_edges is None or not det_edges[slices[:-1]].any(): intens = np.sum(y[slices] * kernel) sigma = np.sqrt(np.sum(esq[slices] * (kernel**2))) - return intens, sigma + status = PEAK_STATUS.STRONG if intens / sigma > weak_peak_threshold else PEAK_STATUS.WEAK + else: + status = PEAK_STATUS.ON_EDGE + intens, sigma = 0.0, 0.0 + return intens, sigma, status def optimise_shoebox(y, esq, current_shape, ipos): @@ -490,7 +600,7 @@ def optimise_shoebox(y, esq, current_shape, ipos): for ncols in lengths[1]: for nbins in lengths[2]: kernel = make_kernel(nrows, ncols, nbins) - intens, sigma = integrate_shoebox_at_pos(y, esq, kernel, ipos) + intens, sigma, _ = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold=0.0) i_over_sig = intens / sigma if i_over_sig > best_i_over_sig: best_i_over_sig = i_over_sig From fb4caaba7931195e0b72eec2079ba2eb6cdd9283 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 15 Nov 2023 16:38:11 +0000 Subject: [PATCH 11/37] Fix bug where expected integer from integer division of a float --- .../PythonInterface/plugins/algorithms/IntegratePeaksSkew.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py index e81109e9106a..43866792cbc5 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksSkew.py @@ -215,7 +215,7 @@ def get_peak_data(self, peak, detid, bank_name, nrows, ncols, nrows_edge, ncols_ """ bank = self.inst.getComponentByName(bank_name) row, col = peak.getRow(), peak.getCol() - drows, dcols = nrows // 2, ncols // 2 + drows, dcols = int(nrows) // 2, int(ncols) // 2 detids, det_edges, irow_peak, icol_peak = self.get_detid_array(bank, detid, row, col, drows, dcols, nrows_edge, ncols_edge) peak_data = PeakData(irow_peak, icol_peak, det_edges, detids, peak, self.ws) return peak_data From 1a04b44a1ac29ff66aacc1cd720b1296c1d70d2c Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 17 Nov 2023 09:33:04 +0000 Subject: [PATCH 12/37] Fix bug trying to integrate peaks on detector edge Note that peaks on the detector edge will not have the shoebox optimised or be used to get the shoebox dimensions for weak peaks --- .../algorithms/IntegratePeaksShoeboxTOF.py | 44 ++++++++++++++++--- 1 file changed, 37 insertions(+), 7 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 256ed43b9c7f..a796d6619228 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -32,6 +32,7 @@ class PEAK_STATUS(Enum): WEAK = "Weak peak" STRONG = "Strong peak" ON_EDGE = "Peak mask is on the detector edge" + ON_EDGE_INTEGRATED = "Peak mask is on the detector edge" NO_PEAK = "No peak detected." @@ -335,7 +336,7 @@ def PyExec(self): nrows = self.getProperty("NRows").value # get these inside loop as overwritten if shoebox optimised ncols = self.getProperty("NCols").value ispec = ws.getIndicesFromDetectorIDs([detid])[0] - itof = ws.binIndexOf(pk_tof, ispec) + itof = ws.yIndexOfX(pk_tof, ispec) # np.clip try/catch bin_width = np.diff(ws.readX(ispec)[itof : itof + 2])[0] # used later to scale intensity if get_nbins_from_b2bexp_params: fwhm = get_fwhm_from_back_to_back_params(peak, ws, detid) @@ -572,14 +573,43 @@ def find_ipks_in_window(ws, peaks, ispecs, ipk, tof_min=None, tof_max=None): def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges=None): - slices = tuple([slice(ii - shape // 2, ii + shape // 2 + 1) for ii, shape in zip(ipos, kernel.shape)]) - if det_edges is None or not det_edges[slices[:-1]].any(): - intens = np.sum(y[slices] * kernel) - sigma = np.sqrt(np.sum(esq[slices] * (kernel**2))) - status = PEAK_STATUS.STRONG if intens / sigma > weak_peak_threshold else PEAK_STATUS.WEAK - else: + slices = tuple( + [ + slice( + np.clip(ii - kernel.shape[idim] // 2, a_min=0, a_max=y.shape[idim]), + np.clip(ii + kernel.shape[idim] // 2 + 1, a_min=0, a_max=y.shape[idim]), + ) + for idim, ii in enumerate(ipos) + ] + ) + + status = PEAK_STATUS.NO_PEAK + 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 + status = PEAK_STATUS.ON_EDGE_INTEGRATED # don't want to optimise weak peak shoeboxes using edge peaks + kernel_slices = [] + for idim in range(len(ipos)): + # start/stop indices in kernel (by default all elements) + istart, iend = 2 * [None] + # 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 + 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)) + kernel = kernel[tuple(kernel_slices)] + # re-normalise the shell to have integral = 0 + inegative = kernel > 0 + kernel[inegative] = -np.sum(kernel[~inegative]) / np.sum(inegative) + intens = np.sum(y[slices] * kernel) + sigma = np.sqrt(np.sum(esq[slices] * (kernel**2))) + if status == PEAK_STATUS.NO_PEAK: + status = PEAK_STATUS.STRONG if intens / sigma > weak_peak_threshold else PEAK_STATUS.WEAK return intens, sigma, status From 2747c0dc5a7575c5ff610daa5d58138fb9fe5dc5 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 17 Nov 2023 11:48:59 +0000 Subject: [PATCH 13/37] Add check for peak TOF in range of workspace x-axis Restructured code to reduce branches a bit --- .../algorithms/IntegratePeaksShoeboxTOF.py | 119 ++++++++++-------- 1 file changed, 66 insertions(+), 53 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index a796d6619228..bfe7e2180751 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -50,7 +50,7 @@ class ShoeboxResult: This class is used to hold data and integration parameters for single-crystal Bragg peaks """ - def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, status, strong_peak=None, ipk_strong=None): + def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, intens_over_sig, status, strong_peak=None, ipk_strong=None): self.peak_shape = list(peak_shape) self.kernel_shape = list(get_kernel_shape(*self.peak_shape)[0]) self.ipos = list(ipos) @@ -67,7 +67,7 @@ def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, status, strong_peak else: self.extents.append([x[0], x[-1]]) # extract peak properties - intens_over_sig = np.round(pk.getIntensityOverSigma(), 1) + intens_over_sig = np.round(intens_over_sig, 1) hkl = np.round(pk.getHKL(), 2) wl = np.round(pk.getWavelength(), 2) tth = np.round(np.degrees(pk.getScattering()), 1) @@ -328,62 +328,70 @@ def PyExec(self): ipks_strong = [] results = np.full(peaks.getNumberPeaks(), None) for ipk, peak in enumerate(peaks): + status = PEAK_STATUS.NO_PEAK + intens, sigma = 0.0, 0.0 + detid = peak.getDetectorID() bank_name = peaks.column("BankName")[ipk] pk_tof = peak.getTOF() - # get shoebox kernel for initial integration - nrows = self.getProperty("NRows").value # get these inside loop as overwritten if shoebox optimised - ncols = self.getProperty("NCols").value - ispec = ws.getIndicesFromDetectorIDs([detid])[0] - itof = ws.yIndexOfX(pk_tof, ispec) # np.clip try/catch - bin_width = np.diff(ws.readX(ispec)[itof : itof + 2])[0] # used later to scale intensity - if get_nbins_from_b2bexp_params: - fwhm = get_fwhm_from_back_to_back_params(peak, ws, detid) - - nbins = max(3, int(nfwhm * fwhm / bin_width)) if fwhm is not None else self.getProperty("NBins").value - nbins = round_up_to_odd_number(nbins) - else: - nbins = self.getProperty("NBins").value - kernel = make_kernel(nrows, ncols, nbins) + # check TOF is in limits of x-axis + xdim = ws.getXDimension() + if xdim.getMinimum() < pk_tof < xdim.getMaximum(): + # get shoebox kernel for initial integration + nrows = self.getProperty("NRows").value # get these inside loop as overwritten if shoebox optimised + ncols = self.getProperty("NCols").value + ispec = ws.getIndicesFromDetectorIDs([detid])[0] + itof = ws.yIndexOfX(pk_tof, ispec) + bin_width = np.diff(ws.readX(ispec)[itof : itof + 2])[0] # used later to scale intensity + if get_nbins_from_b2bexp_params: + fwhm = get_fwhm_from_back_to_back_params(peak, ws, detid) + nbins = max(3, int(nfwhm * fwhm / bin_width)) if fwhm is not None else self.getProperty("NBins").value + nbins = round_up_to_odd_number(nbins) + else: + nbins = self.getProperty("NBins").value + kernel = make_kernel(nrows, ncols, nbins) - # get data array and crop - peak_data = array_converter.get_peak_data( - peak, detid, bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge - ) - x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox) + # get data array and crop + peak_data = array_converter.get_peak_data( + peak, detid, bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge + ) + x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox) - # perform initial integration - intens_over_sig = convolve_shoebox(y, esq, kernel) + # perform initial integration + intens_over_sig = convolve_shoebox(y, esq, kernel) - # identify best shoebox position near peak - ix = np.argmin(abs(x - pk_tof)) - ipos = find_nearest_peak_in_data_window(intens_over_sig, ispecs, x, ws, peaks, ipk, peak_data.irow, peak_data.icol, ix) + # identify best shoebox position near peak + ix = np.argmin(abs(x - pk_tof)) + ipos = find_nearest_peak_in_data_window(intens_over_sig, ispecs, x, ws, peaks, ipk, peak_data.irow, peak_data.icol, ix) - # perform final integration if required - det_edges = peak_data.det_edges if not integrate_on_edge else None - if ipos is None: - status = PEAK_STATUS.NO_PEAK - intens, sigma = 0.0, 0.0 - else: - # integrate at that position (without smoothing I/sigma) - intens, sigma, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges) - if status == PEAK_STATUS.STRONG and do_optimise_shoebox: - kernel, (nrows, ncols, nbins) = optimise_shoebox(y, esq, kernel.shape, ipos) - # re-integrate but this time check for overlap with edge + # perform final integration if required + det_edges = peak_data.det_edges if not integrate_on_edge else None + 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) - - 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) - weak_peaks_list.append(WeakPeak(ipk, ispecs[ipos[0], ipos[1]], x[ipos[-1]], fwhm, ipks_near)) - else: - if status == PEAK_STATUS.STRONG: - ipks_strong.append(ipk) - set_peak_intensity(peak, intens, sigma, bin_width, do_lorz_cor) - if output_file: - # save result for plotting - results[ipk] = ShoeboxResult(ipk, peak, x, y, [nrows, ncols, nbins], ipos, [peak_data.irow, peak_data.icol, ix], status) + if status == PEAK_STATUS.STRONG and do_optimise_shoebox: + kernel, (nrows, ncols, nbins) = optimise_shoebox(y, esq, kernel.shape, ipos) + # 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) + + 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) + weak_peaks_list.append(WeakPeak(ipk, ispecs[ipos[0], ipos[1]], x[ipos[-1]], fwhm, ipks_near)) + 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 + ) + # scale summed intensity by bin width to get integrated area + intens = intens * bin_width + sigma = sigma * bin_width + set_peak_intensity(peak, intens, sigma, do_lorz_cor) if len(ipks_strong): # set function for calculating distance metric between peaks @@ -431,9 +439,13 @@ def PyExec(self): 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) - set_peak_intensity(peak, intens, sigma, bin_width, do_lorz_cor) + # scale summed intensity by bin width to get integrated area + intens = intens * bin_width + sigma = sigma * 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 results[ipk] = ShoeboxResult( ipk, peak, @@ -442,6 +454,7 @@ def PyExec(self): [nrows, ncols, nbins], ipos, [peak_data.irow, peak_data.icol, np.argmin(abs(x - pk_tof))], + i_over_sig, status, strong_peak=strong_pk, ipk_strong=ipk_strong, @@ -639,14 +652,14 @@ def optimise_shoebox(y, esq, current_shape, ipos): return best_kernel, best_peak_shape -def set_peak_intensity(pk, intens, sigma, bin_width, do_lorz_cor): +def set_peak_intensity(pk, intens, sigma, do_lorz_cor): if do_lorz_cor: L = (np.sin(pk.getScattering() / 2) ** 2) / (pk.getWavelength() ** 4) # at updated peak pos else: L = 1 # set peak object intensity - pk.setIntensity(L * intens * bin_width) - pk.setSigmaIntensity(L * sigma * bin_width) + pk.setIntensity(L * intens) + pk.setSigmaIntensity(L * sigma) # register algorithm with mantid From 42bcc4f9789107ae2d20b40af6067a882207ac63 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 Nov 2023 10:33:47 +0000 Subject: [PATCH 14/37] Fix bug for finding nearest peak in window and No Peak plotting --- .../algorithms/IntegratePeaksShoeboxTOF.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index bfe7e2180751..4f1c9e7d2a95 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -53,7 +53,7 @@ class ShoeboxResult: def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, intens_over_sig, status, strong_peak=None, ipk_strong=None): self.peak_shape = list(peak_shape) self.kernel_shape = list(get_kernel_shape(*self.peak_shape)[0]) - self.ipos = list(ipos) + self.ipos = list(ipos) if ipos is not None else ipos self.ipos_pk = list(ipk_pos) self.labels = ["Row", "Col", "TOF"] self.ysum = [] @@ -92,7 +92,8 @@ def plot_integrated_peak(self, fig, axes, norm_func, rect_func): ) # , extent=np.flip(extents, axis=0).flatten()) im.set_norm(norm_func()) # plot shoebox - self.plot_shoebox(iax, ax, rect_func) + if self.ipos is not None: + self.plot_shoebox(iax, ax, rect_func) # plot peak position cen = self.ipos_pk.copy() cen.pop(iax) @@ -552,17 +553,20 @@ def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol tofs = peaks.column("TOF") # get position of nearby peaks in data array in fractional coords shape = np.array(data.shape) - pos_near = [np.array(*np.where(ispecs == ispecs_peaks[ii]), np.array(np.argmin(x - tofs[ii]))) / shape for ii in ipks_near] - # sort data in descending order and select strongest peak nearest to pk (in fractional coordinates) + pos_near = [ + np.r_[np.argwhere(ispecs == ispecs_peaks[ii])[0], np.argmin(abs(x - tofs[ipk_near]))] / shape + for ii, ipk_near in enumerate(ipks_near) + ] + # sort data in descending order and select strongest data point nearest to pk (in fractional coordinates) isort = list(zip(*np.unravel_index(np.argsort(-data, axis=None), data.shape))) pk_pos = np.array([irow, icol, ix]) / np.array(data.shape) imax_nearest = None for ibin in isort: - # calc distance to pk pos + # calc distance to predicted peak position bin_pos = np.array(ibin) / shape - pk_dist_sq = np.sum(pk_pos - bin_pos) + pk_dist_sq = np.sum((pk_pos - bin_pos) ** 2) for pos in pos_near: - if np.sum(pos - bin_pos) > pk_dist_sq: + if np.sum((pos - bin_pos) ** 2) > pk_dist_sq: imax_nearest = ibin break else: From 31a8b524d508b37d172b9486cebf55ef19149d6b Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 Nov 2023 10:45:55 +0000 Subject: [PATCH 15/37] Add threhsold for I/sigma when looking for peak posnearest predicted This provides a shortcut to stop looping over insignificant bins when looking for the peak in a window of data - if no bin can be found with I/sigma > threshold that is nearest to the predicted position than any other peak then it is assinged PEAK_STATUS.NO_PEAK --- .../algorithms/IntegratePeaksShoeboxTOF.py | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 4f1c9e7d2a95..895c4c4fcb43 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -547,7 +547,7 @@ def convolve_shoebox(y, esq, kernel): return intens_over_sig -def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol, ix): +def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol, ix, threshold=2.5): ipks_near, ispecs_peaks = find_ipks_in_window(ws, peaks, ispecs, ipk, tof_min=x.min(), tof_max=x.max()) if len(ipks_near) > 0: tofs = peaks.column("TOF") @@ -562,16 +562,17 @@ def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol pk_pos = np.array([irow, icol, ix]) / np.array(data.shape) imax_nearest = None for ibin in isort: - # calc distance to predicted peak position - bin_pos = np.array(ibin) / shape - pk_dist_sq = np.sum((pk_pos - bin_pos) ** 2) - for pos in pos_near: - if np.sum((pos - bin_pos) ** 2) > pk_dist_sq: - imax_nearest = ibin - break - else: - continue # executed if inner loop did not break (i.e. pixel closer to nearby peak than this peak) - break # execute if inner loop did break and else branch ignored (i.e. found bin closest to this peak) + if data[ibin] > threshold: + # calc distance to predicted peak position + bin_pos = np.array(ibin) / shape + pk_dist_sq = np.sum((pk_pos - bin_pos) ** 2) + for pos in pos_near: + if np.sum((pos - bin_pos) ** 2) > pk_dist_sq: + imax_nearest = ibin + break + else: + continue # executed if inner loop did not break (i.e. pixel closer to nearby peak than this peak) + break # execute if inner loop did break and else branch ignored (i.e. found bin closest to this peak) return imax_nearest # could be None if no peak found else: # no nearby peaks - return position of maximum in data From a537e7af714807f2b5a4f81150079a47cf10ee07 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 Nov 2023 15:15:14 +0000 Subject: [PATCH 16/37] Zero edge of inten_over_sig post convolution where invalid --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 895c4c4fcb43..5c466efcc39f 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -543,11 +543,17 @@ def convolve_shoebox(y, esq, kernel): with np.errstate(divide="ignore", invalid="ignore"): intens_over_sig = yconv / econv intens_over_sig[~np.isfinite(intens_over_sig)] = 0 - intens_over_sig = convolve(intens_over_sig, weights=np.ones(3 * [3]) / (3**3), mode="constant") # 0 pad + ndim = len(y.shape) + intens_over_sig = convolve(intens_over_sig, weights=np.ones(ndim * [ndim]) / (ndim**3), mode="constant") # 0 pad + # zero edges where convolution is invalid + edge_mask = np.ones(intens_over_sig.shape, dtype=bool) + center_slice = tuple(slice(nedge, -nedge) for nedge in (np.asarray(kernel.shape) - 1) // 2) + edge_mask[center_slice] = False + intens_over_sig[edge_mask] = 0 return intens_over_sig -def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol, ix, threshold=2.5): +def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol, ix, threshold=2): ipks_near, ispecs_peaks = find_ipks_in_window(ws, peaks, ispecs, ipk, tof_min=x.min(), tof_max=x.max()) if len(ipks_near) > 0: tofs = peaks.column("TOF") From 52b0a33818835712e55dc30a05986c18c359197e Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 Nov 2023 18:28:38 +0000 Subject: [PATCH 17/37] Improve speed signal.convovle and fewer steps in shoebox optimisation --- .../algorithms/IntegratePeaksShoeboxTOF.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 5c466efcc39f..808917f568b0 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -22,7 +22,8 @@ PropertyCriterion, ) import numpy as np -from scipy.ndimage import convolve +from scipy.ndimage import uniform_filter +from scipy.signal import convolve from IntegratePeaksSkew import InstrumentArrayConverter, get_fwhm_from_back_to_back_params from FindSXPeaksConvolve import make_kernel, get_kernel_shape from enum import Enum @@ -538,13 +539,12 @@ def get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox): def convolve_shoebox(y, esq, kernel): - yconv = convolve(input=y, weights=kernel, mode="nearest") - econv = np.sqrt(convolve(input=esq, weights=kernel**2, mode="nearest")) + yconv = convolve(y, kernel, mode="same") + econv = np.sqrt(convolve(esq, kernel**2, mode="same")) with np.errstate(divide="ignore", invalid="ignore"): intens_over_sig = yconv / econv intens_over_sig[~np.isfinite(intens_over_sig)] = 0 - ndim = len(y.shape) - intens_over_sig = convolve(intens_over_sig, weights=np.ones(ndim * [ndim]) / (ndim**3), mode="constant") # 0 pad + intens_over_sig = uniform_filter(intens_over_sig, size=len(y.shape)) # zero edges where convolution is invalid edge_mask = np.ones(intens_over_sig.shape, dtype=bool) center_slice = tuple(slice(nedge, -nedge) for nedge in (np.asarray(kernel.shape) - 1) // 2) @@ -637,15 +637,15 @@ def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edge return intens, sigma, status -def optimise_shoebox(y, esq, current_shape, ipos): +def optimise_shoebox(y, esq, current_shape, ipos, max_nsteps=25): # find limits of kernel size from data size lengths = [] for idim, length in enumerate(y.shape): max_length = 2 * min(ipos[idim], length - ipos[idim]) + 1 - min_length = max(3, current_shape[idim] // 2) - if not min_length % 2: - min_length += 1 # force odd - lengths.append(np.arange(min_length, max_length - 2 * max(1, max_length // 16), 2)) + min_length = round_up_to_odd_number(max(3, current_shape[idim] // 2)) + step = int(max_length - min_length) // max_nsteps + step = max(2, step + step % 2) + lengths.append(np.arange(min_length, max_length - 2 * max(1, max_length // 16), step)) # loop over all possible kernel sizes best_peak_shape = None best_kernel = None From a136edb7ff2c341fd71911df2439e480ca2d151d Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 23 Nov 2023 17:38:47 +0000 Subject: [PATCH 18/37] Allow center to be optimised as well as shoebox size Use 1D search based on convolution (starting with TOF dimension) --- .../algorithms/IntegratePeaksShoeboxTOF.py | 80 +++++++++++-------- 1 file changed, 48 insertions(+), 32 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 808917f568b0..df7bcbf27974 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -373,7 +373,8 @@ def PyExec(self): # integrate at that position (without smoothing I/sigma) intens, sigma, status = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges) if status == PEAK_STATUS.STRONG and do_optimise_shoebox: - kernel, (nrows, ncols, nbins) = optimise_shoebox(y, esq, kernel.shape, ipos) + 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) @@ -538,18 +539,19 @@ def get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox): return x, y, esq, ispecs -def convolve_shoebox(y, esq, kernel): - yconv = convolve(y, kernel, mode="same") - econv = np.sqrt(convolve(esq, kernel**2, mode="same")) +def convolve_shoebox(y, esq, kernel, mode="same", do_smooth=True): + yconv = convolve(y, kernel, mode=mode) + econv = np.sqrt(convolve(esq, kernel**2, mode=mode)) with np.errstate(divide="ignore", invalid="ignore"): intens_over_sig = yconv / econv intens_over_sig[~np.isfinite(intens_over_sig)] = 0 intens_over_sig = uniform_filter(intens_over_sig, size=len(y.shape)) # zero edges where convolution is invalid - edge_mask = np.ones(intens_over_sig.shape, dtype=bool) - center_slice = tuple(slice(nedge, -nedge) for nedge in (np.asarray(kernel.shape) - 1) // 2) - edge_mask[center_slice] = False - intens_over_sig[edge_mask] = 0 + if mode == "same": + edge_mask = np.ones(intens_over_sig.shape, dtype=bool) + center_slice = tuple(slice(nedge, -nedge) for nedge in (np.asarray(kernel.shape) - 1) // 2) + edge_mask[center_slice] = False + intens_over_sig[edge_mask] = 0 return intens_over_sig @@ -637,30 +639,44 @@ def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edge return intens, sigma, status -def optimise_shoebox(y, esq, current_shape, ipos, max_nsteps=25): - # find limits of kernel size from data size - lengths = [] - for idim, length in enumerate(y.shape): - max_length = 2 * min(ipos[idim], length - ipos[idim]) + 1 - min_length = round_up_to_odd_number(max(3, current_shape[idim] // 2)) - step = int(max_length - min_length) // max_nsteps - step = max(2, step + step % 2) - lengths.append(np.arange(min_length, max_length - 2 * max(1, max_length // 16), step)) - # loop over all possible kernel sizes - best_peak_shape = None - best_kernel = None - best_i_over_sig = -np.inf # null value - for nrows in lengths[0]: - for ncols in lengths[1]: - for nbins in lengths[2]: - kernel = make_kernel(nrows, ncols, nbins) - intens, sigma, _ = integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold=0.0) - i_over_sig = intens / sigma - if i_over_sig > best_i_over_sig: - best_i_over_sig = i_over_sig - best_peak_shape = (nrows, ncols, nbins) - best_kernel = kernel - return best_kernel, best_peak_shape +def optimise_shoebox(y, esq, peak_shape, ipos, nfail_max=2): + best_peak_shape = list(peak_shape) + best_ipos = list(ipos) + for idim in range(3)[::-1]: + nfailed = 0 + max_intens_over_sig = -np.inf + this_peak_shape = best_peak_shape.copy() + this_peak_shape[idim] = min(1, round_up_to_odd_number(peak_shape[idim] // 2) - 2) + while True: + # make kernel + this_peak_shape[idim] = this_peak_shape[idim] + 2 + kernel = make_kernel(*this_peak_shape) + # get slice the size of kernel along all other dims centered on best ipos from other dims + slices = [ + slice(np.clip(ii - kernel.shape[idim] // 2, 0, y.shape[idim]), np.clip(ii + kernel.shape[idim] // 2 + 1, 0, y.shape[idim])) + for idim, ii in enumerate(best_ipos) + ] + # incl. elements along dim that give valid convolution in initial peak region + slices[idim] = slice( + np.clip(ipos[idim] - peak_shape[idim] // 2 - kernel.shape[idim] // 2, 0, y.shape[idim]), + np.clip(ipos[idim] + peak_shape[idim] // 2 + kernel.shape[idim] // 2 + 1, 0, y.shape[idim]), + ) + slices = tuple(slices) + if y[slices].shape[idim] <= kernel.shape[idim]: + # no valid convolution region + break + this_intens_over_sig = np.squeeze(convolve_shoebox(y[slices], esq[slices], kernel, mode="valid")) + imax = np.argmax(this_intens_over_sig) + if this_intens_over_sig[imax] > max_intens_over_sig: + max_intens_over_sig = this_intens_over_sig[imax] + best_ipos[idim] = imax + max(ipos[idim] - peak_shape[idim] // 2, kernel.shape[idim] // 2) + best_peak_shape[idim] = this_peak_shape[idim] + nfailed = 0 + else: + nfailed += 1 + if not nfailed < nfail_max: + break + return best_ipos, best_peak_shape def set_peak_intensity(pk, intens, sigma, do_lorz_cor): From 03f15b3866f5511d40434504e022ffd231b15232 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 24 Nov 2023 10:38:40 +0000 Subject: [PATCH 19/37] Fix TOF axes limits in output pdf plots --- .../algorithms/IntegratePeaksShoeboxTOF.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index df7bcbf27974..87baab262016 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -58,15 +58,11 @@ def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, intens_over_sig, st self.ipos_pk = list(ipk_pos) self.labels = ["Row", "Col", "TOF"] self.ysum = [] - self.extents = [] + self.xmin, self.xmax = np.round(x[0], 1), np.round(x[-1], 1) self.status = status # integrate y over each dim for idim in range(len(peak_shape)): self.ysum.append(y.sum(axis=idim)) - if idim < 2: - self.extents.append([0, y.shape[idim]]) - else: - self.extents.append([x[0], x[-1]]) # extract peak properties intens_over_sig = np.round(intens_over_sig, 1) hkl = np.round(pk.getHKL(), 2) @@ -86,11 +82,7 @@ def __init__(self, ipk, pk, x, y, peak_shape, ipos, ipk_pos, intens_over_sig, st def plot_integrated_peak(self, fig, axes, norm_func, rect_func): for iax, ax in enumerate(axes): - extents = self.extents.copy() - extents.pop(iax) - im = ax.imshow( - self.ysum[iax], aspect=self.ysum[iax].shape[1] / self.ysum[iax].shape[0] - ) # , extent=np.flip(extents, axis=0).flatten()) + im = ax.imshow(self.ysum[iax], aspect=self.ysum[iax].shape[1] / self.ysum[iax].shape[0]) im.set_norm(norm_func()) # plot shoebox if self.ipos is not None: @@ -104,6 +96,14 @@ def plot_integrated_peak(self, fig, axes, norm_func, rect_func): labels.pop(iax) ax.set_xlabel(labels[1]) ax.set_ylabel(labels[0]) + # overwrite axes tick labels for TOF axes + nticks = 4 + xticklabels = nticks * [""] + xticklabels[0] = str(self.xmin) + xticklabels[-1] = str(self.xmax) + for ax in axes[:-1]: + ax.set_xticks(np.linspace(*ax.get_xlim(), nticks)) + ax.set_xticklabels(xticklabels) # set title fig.suptitle(self.title) From 747a73d1ad382769fa8cc395160a1bd414054ef8 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 24 Nov 2023 10:53:34 +0000 Subject: [PATCH 20/37] Move sqrt of error conv to with statement to supress warnings Sometimes 0s can be infinitessimally negative which causes this warning --- .../algorithms/IntegratePeaksShoeboxTOF.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 87baab262016..e434ff07ab75 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -449,18 +449,10 @@ def PyExec(self): 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( - ipk, - peak, - x, - y, - [nrows, ncols, nbins], - ipos, - [peak_data.irow, peak_data.icol, np.argmin(abs(x - pk_tof))], - i_over_sig, - status, - strong_peak=strong_pk, - ipk_strong=ipk_strong, + ipk, peak, x, y, peak_shape, ipos, ipos_predicted, i_over_sig, status, strong_peak=strong_pk, ipk_strong=ipk_strong ) elif weak_peak_strategy == "NearestStrongPeak": raise ValueError( @@ -539,10 +531,10 @@ def get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox): return x, y, esq, ispecs -def convolve_shoebox(y, esq, kernel, mode="same", do_smooth=True): +def convolve_shoebox(y, esq, kernel, mode="same"): yconv = convolve(y, kernel, mode=mode) - econv = np.sqrt(convolve(esq, kernel**2, mode=mode)) with np.errstate(divide="ignore", invalid="ignore"): + econv = np.sqrt(convolve(esq, kernel**2, mode=mode)) intens_over_sig = yconv / econv intens_over_sig[~np.isfinite(intens_over_sig)] = 0 intens_over_sig = uniform_filter(intens_over_sig, size=len(y.shape)) From 9a7c27bbd37f8763b2e79fe9dd417f65e4f2981c Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 24 Nov 2023 11:08:30 +0000 Subject: [PATCH 21/37] Only try plotting if output file specified --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index e434ff07ab75..12a5b807b776 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -462,7 +462,8 @@ def PyExec(self): ) # plot output - plot_integration_reuslts(output_file, results) + if output_file: + plot_integration_reuslts(output_file, results) # assign output self.setProperty("OutputWorkspace", peaks) From ab07a5f34bbec420e373d2b792d997f9198ebbca Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 24 Nov 2023 11:22:19 +0000 Subject: [PATCH 22/37] Add progress reporter --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 12a5b807b776..ac2e09707ef1 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -12,6 +12,7 @@ WorkspaceUnitValidator, FileProperty, FileAction, + Progress, ) from mantid.kernel import ( Direction, @@ -329,7 +330,9 @@ def PyExec(self): weak_peaks_list = [] ipks_strong = [] results = np.full(peaks.getNumberPeaks(), None) + prog_reporter = Progress(self, start=0.0, end=1.0, nreports=peaks.getNumberPeaks()) for ipk, peak in enumerate(peaks): + prog_reporter.report("Integrating") status = PEAK_STATUS.NO_PEAK intens, sigma = 0.0, 0.0 @@ -463,7 +466,8 @@ def PyExec(self): # plot output if output_file: - plot_integration_reuslts(output_file, results) + prog_reporter.resetNumSteps(int(len(results) - np.sum(results is None)), start=0.0, end=1.0) + plot_integration_reuslts(output_file, results, prog_reporter) # assign output self.setProperty("OutputWorkspace", peaks) @@ -486,7 +490,7 @@ def round_up_to_odd_number(number): return number -def plot_integration_reuslts(output_file, results): +def plot_integration_reuslts(output_file, results, prog_reporter): # import inside this function as not allowed to import at point algorithms are registered from matplotlib.pyplot import subplots, close from matplotlib.patches import Rectangle @@ -497,6 +501,7 @@ def plot_integration_reuslts(output_file, results): with PdfPages(output_file) as pdf: for result in results: if result: + prog_reporter.report("Plotting") fig, axes = subplots(1, 3, figsize=(12, 5), subplot_kw={"projection": "mantid"}) fig.subplots_adjust(wspace=0.3) # ensure plenty space between subplots (want to avoid slow tight_layout) result.plot_integrated_peak(fig, axes, LogNorm, Rectangle) From 0246b7f3db813907c928c51fdb9f3a90fcd8883c Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 24 Nov 2023 11:42:35 +0000 Subject: [PATCH 23/37] Add unit test --- .../python/plugins/algorithms/CMakeLists.txt | 1 + .../IntegratePeaksShoeboxTOFTest.py | 192 ++++++++++++++++++ 2 files changed, 193 insertions(+) create mode 100644 Framework/PythonInterface/test/python/plugins/algorithms/IntegratePeaksShoeboxTOFTest.py diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt index da51be20860b..7b16f9e4e25c 100644 --- a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt +++ b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt @@ -69,6 +69,7 @@ set(TEST_PY_FILES IndexPeaksTest.py IndirectTransmissionTest.py IndexSatellitePeaksTest.py + IntegratePeaksShoeboxTOFTest.py IntegratePeaksSkewTest.py LeadPressureCalcTest.py LoadAndMergeTest.py diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/IntegratePeaksShoeboxTOFTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/IntegratePeaksShoeboxTOFTest.py new file mode 100644 index 000000000000..7a29b1026ea1 --- /dev/null +++ b/Framework/PythonInterface/test/python/plugins/algorithms/IntegratePeaksShoeboxTOFTest.py @@ -0,0 +1,192 @@ +# 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 mantid.simpleapi import ( + IntegratePeaksShoeboxTOF, + CreatePeaksWorkspace, + AddPeak, + LoadParameterFile, + AnalysisDataService, + SortPeaksWorkspace, +) +from testhelpers import WorkspaceCreationHelper +from numpy import array, sqrt +import tempfile +import shutil +from os import path + + +XML_PARAMS = """ + + + + + + + + + + + + + + +""" + + +class FindSXPeaksConvolveTest(unittest.TestCase): + @classmethod + def setUpClass(cls): + # load empty instrument with RectangularDetector banks and create a peak table + cls.ws = WorkspaceCreationHelper.create2DWorkspaceWithRectangularInstrument(1, 7, 13) # nbanks, npix, nbins + AnalysisDataService.addOrReplace("ws_rect", cls.ws) + axis = cls.ws.getAxis(0) + axis.setUnit("TOF") + # fake peak centred on ispec=30 (detid=79) and TOF=5 - near middle of bank + peak_1D = array([0, 0, 4, 6, 4, 0, 0, 0, 0, 0, 0, 0, 0]) + cls.ws.setY(30, cls.ws.readY(30) + peak_1D) + for ispec in [23, 29, 30, 31, 37]: + cls.ws.setY(ispec, cls.ws.readY(ispec) + peak_1D) + # fake peak centred on ispec=12 (detid=61) and TOF=7 - near detector edge + cls.ws.setY(12, cls.ws.readY(12) + peak_1D[::-1]) + for ispec in [5, 11, 12, 13, 19]: + cls.ws.setY(ispec, cls.ws.readY(ispec) + peak_1D[::-1]) + # add background and set errors + for ispec in range(cls.ws.getNumberHistograms()): + cls.ws.setY(ispec, cls.ws.readY(ispec) + 0.5) + cls.ws.setE(ispec, sqrt(cls.ws.readY(ispec))) + # add peaks + cls.peaks = CreatePeaksWorkspace(InstrumentWorkspace=cls.ws, NumberOfPeaks=0, OutputWorkspace="peaks") + AddPeak(PeaksWorkspace=cls.peaks, RunWorkspace=cls.ws, TOF=5, DetectorID=85) + AddPeak(PeaksWorkspace=cls.peaks, RunWorkspace=cls.ws, TOF=9, DetectorID=69) + [cls.peaks.getPeak(ipk).setHKL(ipk, ipk, ipk) for ipk in range(cls.peaks.getNumberPeaks())] + # Add back-to-back exponential params + LoadParameterFile(cls.ws, ParameterXML=XML_PARAMS) + # output file dir + cls._test_dir = tempfile.mkdtemp() + + @classmethod + def tearDownClass(cls): + AnalysisDataService.clear() + shutil.rmtree(cls._test_dir) + + def _assert_found_correct_peaks(self, peak_ws, i_over_sigs=[6.4407, 2.4503]): + self.assertEqual(peak_ws.getNumberPeaks(), 2) + peak_ws = SortPeaksWorkspace( + InputWorkspace=peak_ws, OutputWorkspace=peak_ws.name(), ColumnNameToSortBy="DetID", SortAscending=False + ) + pk = peak_ws.getPeak(0) + self.assertEqual(pk.getDetectorID(), 85) + self.assertAlmostEqual(pk.getTOF(), 5.0, delta=1e-8) + self.assertAlmostEqual(pk.getIntensityOverSigma(), i_over_sigs[0], delta=1e-4) + pk = peak_ws.getPeak(1) + self.assertEqual(pk.getDetectorID(), 69) + self.assertAlmostEqual(pk.getTOF(), 9.0, delta=1e-8) + self.assertAlmostEqual(pk.getIntensityOverSigma(), i_over_sigs[1], delta=1e-4) + + def test_exec_IntegrateIfOnEdge_False(self): + out = IntegratePeaksShoeboxTOF( + InputWorkspace=self.ws, + PeaksWorkspace=self.peaks, + OutputWorkspace="peaks1", + GetNBinsFromBackToBackParams=False, + NRows=3, + NCols=3, + NBins=3, + WeakPeakThreshold=0.0, + OptimiseShoebox=False, + IntegrateIfOnEdge=False, + ) + # check edge peak integrated but lower I/sigma as shell would overlap edge + self._assert_found_correct_peaks(out, 2 * [0.0]) + + def test_exec_IntegrateIfOnEdge_True(self): + out = IntegratePeaksShoeboxTOF( + InputWorkspace=self.ws, + PeaksWorkspace=self.peaks, + OutputWorkspace="peaks2", + GetNBinsFromBackToBackParams=False, + NRows=3, + NCols=3, + NBins=3, + WeakPeakThreshold=0.0, + OptimiseShoebox=False, + IntegrateIfOnEdge=True, + ) + # check edge peaks integrated but lower I/sigma on second peak as kernel overlaps detector edge + self._assert_found_correct_peaks(out) + + def test_exec_GetNBinsFromBackToBackParams(self): + out = IntegratePeaksShoeboxTOF( + InputWorkspace=self.ws, + PeaksWorkspace=self.peaks, + OutputWorkspace="peaks3", + GetNBinsFromBackToBackParams=True, + NRows=3, + NCols=3, + WeakPeakThreshold=0.0, + OptimiseShoebox=False, + IntegrateIfOnEdge=True, + ) + # check edge peaks integrated but lower I/sigma on second peak as kernel overlaps detector edge + self._assert_found_correct_peaks(out) + + def test_exec_OptimiseShoebox(self): + # make kernel larger than optimum + out = IntegratePeaksShoeboxTOF( + InputWorkspace=self.ws, + PeaksWorkspace=self.peaks, + OutputWorkspace="peaks4", + GetNBinsFromBackToBackParams=False, + NRows=5, + NCols=5, + NBins=3, + WeakPeakThreshold=0.0, + OptimiseShoebox=True, + IntegrateIfOnEdge=True, + ) + # check I/sigma first peak unchanged, second peak improved! + self._assert_found_correct_peaks(out, i_over_sigs=[6.4407, 4.0207]) + + def test_exec_OptimiseShoebox_respects_WeakPeakThreshold(self): + # make kernel larger than optimum + out = IntegratePeaksShoeboxTOF( + InputWorkspace=self.ws, + PeaksWorkspace=self.peaks, + OutputWorkspace="peaks5", + GetNBinsFromBackToBackParams=False, + NRows=5, + NCols=5, + NBins=3, + WeakPeakThreshold=10.0, + OptimiseShoebox=True, + IntegrateIfOnEdge=True, + ) + # check I/sigmas much worse if not optimised + self._assert_found_correct_peaks(out, i_over_sigs=[2.7001, 1.1531]) + + def test_exec_OutputFile(self): + out_file = path.join(self._test_dir, "out.pdf") + IntegratePeaksShoeboxTOF( + InputWorkspace=self.ws, + PeaksWorkspace=self.peaks, + OutputWorkspace="peaks1", + GetNBinsFromBackToBackParams=False, + NRows=3, + NCols=3, + NBins=3, + WeakPeakThreshold=0.0, + OptimiseShoebox=False, + IntegrateIfOnEdge=False, + OutputFile=out_file, + ) + # check output file saved + self.assertTrue(path.exists(out_file)) + + +if __name__ == "__main__": + unittest.main() From dff93c29dc34c0493d1a484674dab79080d92a21 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 28 Nov 2023 13:49:06 +0000 Subject: [PATCH 24/37] Add algorithm doc page --- .../IntegratePeaksShoeboxTOF-v1.rst | 89 ++++++++++++++++++ .../IntegratePeaksShoeboxTOF_OutputFile.png | Bin 0 -> 66956 bytes 2 files changed, 89 insertions(+) create mode 100644 docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst create mode 100644 docs/source/images/IntegratePeaksShoeboxTOF_OutputFile.png diff --git a/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst new file mode 100644 index 000000000000..ca32e8a34ee0 --- /dev/null +++ b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst @@ -0,0 +1,89 @@ +.. algorithm:: + +.. summary:: + +.. relatedalgorithms:: + +.. properties:: + +Description +----------- + +This is an algorithm to integrate single-crystal Bragg peaks in a :ref:`MatrixWorkspace ` by +summing up counts in a shoebox of size ``NRows`` and ``NCols`` on the detector and ``NBins`` along the Time-of-flight +(TOF) direction which can be specified in one of two ways: + +1. Provide ``NBins`` - number of TOF bins in the kernel + +2. Setting ``GetNBinsFromBackToBackParams=True`` and providing ``NFWHM`` - in which case ``NBins`` will be NFWHM x FWHM + of a :ref:`BackToBackExponential ` peak at the center of each detector panel/bank at the + middle of the spectrum. + +Note to use method 2, back-to-back exponential coefficients must be defined in the Parameters.xml file for the +instrument. + +The integration requires a background shell with negative weights, there are approximately the same number of bins in +the background shell as in the peak region. + +The algorithm proceeds as follows: + +1. Take a window of the data, ``NShoeboxInWindow`` the size of the shoebox (``NRows`` x ``NCols`` x ``NBins``) + +2. Convolve the shoebox kernel with the data in the window to find the position with largest Intensity/sigma + (closest to the predicted peak position than to any other peaks in the table). + +3. Integrate using the shoebox at the optimal position + +4. If the peak is strong (Intensity/sigma < ``WeakPeakThreshold``) and ``OptimiseShoebox=True`` then optimise the shoebox + dimensions to maximise Intensity/sigma + +5. If the peak is weak, it can be integrated using the initial shoebox (``WeakPeakStrategy="Fix"``) or using the + shoebox dimensions from the nearest strong peak (``WeakPeakStrategy="NearestStrongPeak"``) + +When looking for the nearest strong peaks for ``WeakPeakStrategy="NearestStrongPeak"``, the algorithm first checks for +peaks in detector IDs in the data window around the peak, if one isn't found the nearest peak search depends on the +parameter ``GetNBinsFromBackToBackParams``. + +1. If ``GetNBinsFromBackToBackParams=True`` then the closest peak is defined as the one with the smallest angle between + QLab vectors and the TOF extent of the shoebox is scaled by the ratio of the FWHM of the weak and strong peak. + +2. If ``GetNBinsFromBackToBackParams=False`` then then the closest peak is defined as the one nearest in QLab. + +Optionally if ``OutputFile`` is provided a pdf can be output that shows the shoebox kernel and the data integrated along +each dimension like so + +.. figure:: ../images/IntegratePeaksShoeboxTOF_OutputFile.png + :align: center + :width: 50% + :alt: (Left) Found fractional TOF window and estimated curves at 4 different wavelengths from linear fit shown on + (Right) + + +Usage +----- + +**Example - IntegratePeaksShoeboxTOF** + +.. testcode:: exampleIntegratePeaksShoeboxTOF + + from mantid.simpleapi import * + + Load(Filename="SXD23767.raw", OutputWorkspace="SXD23767") + CreatePeaksWorkspace(InstrumentWorkspace="SXD23767", NumberOfPeaks=0, OutputWorkspace="peaks") + AddPeak(PeaksWorkspace="peaks", RunWorkspace="SXD23767", TOF=8303.3735339704781, DetectorID=7646) + + peaks_out = IntegratePeaksShoeboxTOF(InputWorkspace="SXD23767", PeaksWorkspace="peaks", + GetNBinsFromBackToBackParams=True, WeakPeakThreshold=0.0, LorentzCorrection=False) + + print(f"I/sigma = {peaks_out.getPeak(0).getIntensityOverSigma():.2f}") + +**Output:** + +.. testoutput:: exampleIntegratePeaksShoeboxTOF + + I/sigma = 100.49 + + +.. categories:: + +.. sourcelink:: diff --git a/docs/source/images/IntegratePeaksShoeboxTOF_OutputFile.png b/docs/source/images/IntegratePeaksShoeboxTOF_OutputFile.png new file mode 100644 index 0000000000000000000000000000000000000000..51b7a5997ef7fd24ababe99407c400248b310e4d GIT binary patch literal 66956 zcmeFZ2UJwq)-KwuZE8S8a&9FE2og#tk<_hFh=52C2@;DckW?Z;GOZ#=EXfooC`ob- zf&>Af2n8iMgF=&}$WUasYXf#T_w+e^`u*>Y``$QXw5=l4UTdy7zd7eOzqNgD-Bh6Z zp85N4zWIhq@w%M)H{X00`pq}rE*?J$zG66v9RR<5Yp1So^_#2~)+z9p@63Ku`RSW) za}j^2 zYZ$-1>V8Q6gr8fuOXBm8jUQudZftC{W)=i$3*(Q{YRU=K`)H|gyh(eNj94&Fr!8=w znM`Jr%(52lqj^&pW+SbByN)nWmpfB&_d{2Svb`K0WyN`gOft(1pWY(p#}QUc$KMSV zl81H;Kgg0L*2hCJ6}WZq-|(F7CW_(e@v~m1_y5NgjN|(H=i_JF>Xpy@_K){Wc~1QL zjn6KLhC<(HOE2ZIwKmdg7FX^xjx)Q+F70B!(556b6STav+^%XUhutiRJqr{XtQL=0iYWWH||S$9`_8+-zW7YB= zgZEZGg_h{@d~bcvLQOE#bk{!YSo55w3Ngs9F%js(3%bWXi?I(nV2V?mOqdgU*pj1#;fdr#6ScmOT%1Nhl;TZG@F9MF8V%Kb@c72zOpXBvy77z`MkWvkd zaXRTvP=P+0(|qg)=t-f*${z1=+--}jJVB{qT8 zzhcxp(5^>Mxe&Rxcd-L2|I(u;Uc1aIicIYOMVO|+#FAs>i?nJ(P4Sa}Mf&IE(!Wee zC`imU!Gqj`79&DmDnQSRKX!1U6r@0Am=xi(A0b=-off5d#yz7cN*c=~uH`Z9>TOdk z;dd6bu=hCWgLBHPCQ!o^2?7#F$NaMtC|?NRWfYt*Gz7pQBR?}ZSX7rWBBZ=?v{Kgm z1S#|N8yTS7!(NyIW1ssuZwCv%{61 z+;yg7a^5#FJM+TS9V2bHK|tu`U)%bJxZraN^u8 zJJ#bgO$!*-m2r&PNJhY71kJGH9?UTQw7ccSn(EUY+>1NMoA;W#mjRlU)F~%17&sX(MYQz~?~vU$6l6DARZ7*`@)lQq7!y2)aSM)? z(BhroJF1G=xfZr%Kz0mO=7!c0a1aV+k8mutq)(B`F|U#L$kdk?mt=RBHn-3}IL9No z1=Axww}G`ahw)GF;a$OQlrL}Azrtx~san>U5LP^J#n1uEeKqmp#q#Srl8*?Kx%Ru= z_(nZNGD6+UT{VOa8<5OB=?t6FM5}w%Mm*4U*h>#7XZKyr?5CNg;^w<~_Pg{dBK{<# z6?L}f_B0%iMe(eZ;lj3ibTK<>jfUR&OM6JL{B`L+nI?u|E1aJ+@YdZ}m*J|$zF)|G z8(xs1i}^-JT2OU*CU#Nlc9;+5IGxK9Ry{x%7QS~BE?<56m0RQ~suznzA-m6F)B?im zmH9ayxlTw9mnij`c;Cj7gi<=|%W#^MZ|V#9Jio-@p-&ki=ULK1{xK=P?^icZ?vkWs zTeOt1^&Na1$jt?g-39$d`b}*;zG-a^TbSPhO}V`0CkLE}R#tNbM_QCwz+AJT@U^D4 z-kbC^?o~w<`$owB3e0>*k5#($2b@+)FZV?nmgy@87gui2Fvu>@$l5Lx6_;Dbp)uzL z199SGnMo%0y~Ggu@+8e6BYx*5h2F;9^3Uxt!@T%){To^wO8IeUW}Om&ram!#B`B6g zP9H?#1HVTmbqZCK&_kq=!%t0l)I}8hNP*n2h@eMI!aHu7&1?J|-3bK7H1}Xf%udMp zjb}z&OnYf5>g#ETm<+IPG@;09;+rWk4Li+312X3cWg_`_84pDJ8w6(OnYwr3Yxjwe z?Eo23UdZ&VPxf)=1vi?jxBv29n2nU}JNG^o{gvtfa$n$d%=@@-{dMIbo)4vyu{0o1 zvi6NAN3$8aNWRla%1)Pdn*WjI?%AM#SqboS6G8h6(9Bh_>iyVvbaklQl>y4yY;;!z z0{6ZJ{_XzfrYw^3i8kUzLgmYC1p(# zRav#-MUEbu#qP`FJKLM@5|*^Ld$b&bC^^>a6~x_jmubMkzvG-X+9vh)Ua)wm#8EOt z?@`vgyeg7wS)5z&%Nn16!|)1Ob!z&0!nV|`;dmAl|1>VRpY}MiygyQWz z8Zz;k)O@Ia0G+FnSSyF>#T!bgRm8uE^Kgl@r1Td;>AbWWl~Gz>J2|am8Z*&{*kbVb z%pg9lM!Wj?l&4eV^jZ%>Htle&BS-lVz4(v?s;Dtv{L)f$U#oISVa^quK{P(W%y-kK zz%zbJ9@gy~;Cd6Bs`i;ZgU)^@R5o8-J8ihT(L7GlPQKPV(2rMdjG;zv^`Hx^Cxj?WZo zKE1vkx4Z5S#mMI6&0S_~S%kH!LilyzBUf8?aM9*9i$+S4>7Y&HzzVnuMV2;h9Q+6WE@zl`6lr0SxWMt;+wt~ z;eoHqbMHSh;gQGO*LJffF?QF&6W@wp#Rt)}9pb-~(!}_=M|GJ0+!P|8?uthUZbpyk$sI*)-gJyyVBMytzS5l;DYkR zy?wtpI4R1l*~CJblc@Vo%wv|Wtx<9F#UOX7qeF%Rj`lD2jZ52Q(T@qHENUr86ypxD z$365mxjHX}cG8Z|5p|m6wBg*= zt@BM~6iE)L=!#J+U)2xP+6`NbES+=N3&4@Ey&$$d&x@7|E>HU9^Fn|Hi&X!Ih(K<5 zmo$eWKZgUKt(E~_lZQ+Dt6x^~nz!o5i^{f57M?D?RMBEkP1Uo{IA#gEfES-HL%90a zwKZSLt61Xau(RC04VN5B7@{g;y3Q@WS$Av3MUS@T(eP{eK)2;g1g0`VJZG9u#8ID> zsKvuKrNg#VvBdf%?!n3I;LWU62w%?LJn)P|9*X&VE>E7~VQJR%*7?&*eaSft{XcFz zunJ{mqAThx1WwUCVs+U19z9Yve;o0;OsyJj)mV4MB|otYC$gmFT^li7U>CeN-g*|7 ztI63otUgy57Utt!k zqOEt;0IUOdo;m*b)Dkd?4TLl6r7wl=J>qSt8RLX3I;@qTodco#Xp zqp4%5uJ-I(+w9K#7x4J4@Z!SFFk09T`BTKZ7>Z~v_Gt!5be@b}%!AqE?`$N+A2cT{ z4?KL@yCdQ<`%YuP4xG~m&~5+}9I$A_Cy-|7Y6wYr5=&vr87|_}(wKlKHB;vVtWvK; z(>g4ssgLE9n4gZrbtNXcf$RNsvZXv8!?(T4QD}TT_F{-SLQA?^TiJhdRG-VP+=tg) z7sikHQb{{DKwZnYfXw}1>X5P7q8;xU5W8tJvme5T_DvwT<-{VpqbY8-v*}=EHI zZ3$~ep*1W~dPAVrz&wm14g2u&_=GzjkFZ`T41Bw*)0n?YG+|j;JnFM;Wj&x&E_LbH zsqLpSWPL>M4(eT@iAz@TQzk7f@0F+cMgfP00BjN;B7jLD_Q5Ib{xid~KK2=#HmN%z z;57PHUw|x`OS^A}8w9BleE)T4y#^WfLaw4mLE^6>w8VQo>~zEf27ZWy(+Si?-Wk#r zrL&kQgfrbRki?&DL|hTg{a$OB7t6%z7)0u)(HhQEEQ(XxR-95I<8TT)&G{O)HWE>+ z%LduhfmAgnnXI9nG=7*8fuBGmE$QjaGYeJ>92RehF!kMRStOH*XI2LZzsDIH?y&^b z;D%cjInMVESBqz8Wo;LzBw8P`KCCTXz_VL)J*=fUrB)a0x0e)6+I>1%cIl-kMklOT z|J-5!u)s|CcQ@EhTE|jUc+k()7VDAxixf>F;IAg}NTfRDYeOEhrtT=A03;;9f z56fj6Q!3iV_kTY&iXm$?g@maXa`g)%(AOkpvAuQSF~*H8*Q`ca$6#!IADJD?7^!`f z4>Rj0c8FP;gatY&mED|Hm{yI$Z_H~VLX2lnD>>7cR8?si&l=k++P~L41 zbGHED_G~SzsN>88vw-t3N@C0%DjyCqfmc+#+g-eP=Lq?g2OU&%QSH}8Q6`ziKYl1p zF@(AkwQK*jef;aSovpor%T!WxT(l>-*GwqZc^ubF#A7)yZMDwxfYHz|BcX)eIc-e6 zdB_LW<8@l(prF+x$ccu=^GMz#5nE$kFU)_g(f89nyuaG7UB5!jkPLrvFpd_@io0hl z%OT}--xjb#%`rFRlI!l>dI$Aq=lu5zL{=&o58;*5m((|oe0OrOCF|JIhK|cr77xMq zTs1vtAQIFFT$@3Ol;pZLh=UR#cZF>CH|E5V7vQ1W{X^(H&;EXuh;B_`yT%r?qgoPZ zk)2U{M7orJ@g;0HhNCKfBRp@$Eu#BN^tKZ_(He7qVL{9@Gb&Ng6?W*=avm?Km*-Xxk)MLw{JCC11N-OAgD9F?hj8RxCeSILSj zR~fqHN8mEZ7BYm)`|yA^xXesOD~kvbSaXZ-lc-+z>A?}E7PMtQsZ$1P zWuW>&DPL;rMPkRBf_GzrD>-ojdwj`5mi?2v%wJnJFZ!gKa1sj_yDotqS2VQ=D>gSz zqixiaCI{ynYFeUQ@?PN`Q{cc|;Lz~@!lu?@z1VHCD}-;1Q3H4qw<1JK)661AEfAB; z-apu)Lzy<|206~#{)2X_l@HSpC89Ru$d~56oPpx$q>NEVFOO^ZCVdW(U>H8hKc;}t za+XY8ZBFbf)GlQ!c0))M7rw@!94i*{8FhHnhvf3E=s1s=HRv^p-*E4dt-t^eIewOM ze+SxZuG4SI#*qA&VMf8wvFG1vH2MoaY2d(w-&tIW!k$;#q4&wEZJTj&Zc2Nx`TTpW zxJAd{R7S<*kelt99dGIv<&$yH$9e6$<%{p4CDEp|x;#y14=0B);C{HI^gYgQK_%e~ zN4TMkMctO*NqzsM{)I(0)v*Dxr+;SZts$+v2X5~E!T2pQMkya2>ep+!8}WLNkHvod zP*;`)N(v1TpF--5$zzzqu|La?Mb4R~bv9pH;)%cOKL*pC z2+Ao7IG+g*&DcTzKv&b?(Rd4DTfX}m8OFmCG&m-2ceiv#_TW@wZyz6T1vwFW!Z)84 zCP>4z0l8P`bjQMxZtc&8{WY_i0_b5*0)rD6wMKY2Hjl1H{Aoy+=);~M?x7e}fp>SG zeD1GqQR^cdJmo9j5ai87QyUv__FrB$U=5qZ%LDc0Jboxv(+OygI((3 zdQ?mjOqH)+(-oHb-mXuz=p}RtvlkEG<`RMOwZsStp=x?d`vqKvK)P(o!-Re=H&JjY z`aB}N>;XZU*2Ft}%&Y|65mp@Pw!`8YxLNVWtznbe!xes{JAF3ay28 z)dSKS$~{`9Q2VtGm_{|bi%%eb_? zyJOjPC{j&M91!XDrNwwXRVSS3iUpQRr|9IJ^JV zx2gM2^_qN7$+hc#lf2uvD%*KLImG}pfH3dbPoz=ngyC&dTl(RR3eLSR*~$G6B0uzr zzjgI{KL4Ndl#nIBr~HKx!+#(5f797@d_oDi;VZpagKZI*) zK+F#vakRpCc&lxig8Vv#SAfv^8mtgcM5Yv;2?`c89wmWCp_C8xr=t49qDsABQ6!BlQ%hIRgmW3QZ+RMRT>$kH@J_{=m zPIa|lwyR}j1ju74@elGhRI)e7eDWyrqzc9f9~OUQg*UA|I{K^ymR)`(H#t6cAVDWy zw7V_FU_sD0Zkd-^Ri~M7%1koRq#nS#t}{+y0a~hAh1k$cs@6>?#bl(i6enu@o?!AM5<&Y|8VJfovs2ba1Vcb(%NY z-38_0+=E)ehpwC+gRRRMzN62DiIa8-M~Qr}sSW1Gj#FDd25~@LAUupQ>{ctK4yy@90qS6*D7zC2G~JnDpNdpQmupC zlApOrI+~l`POw2)_nF;`)7LxJgDp&#sBwx5M<4w)cT+&OIYZTHb78MtIaJYFw~?HO zO?Zmxu1jv(YSZA)8ppOyAgE$B&}}PAjQh@I%F}7Atd+1We2AzNH7j;#RAF{E*sM#% zh`e60yk%tDiWa=E?@NB->(YLWB((y?SvAc!?D7t}8Zfz;)`S%4F3n33u9uWF$8|$+ z1*a!2pM9RP#0dqHXy^#lp`}1!mN;vgI3e{S4yi&n<<&ImF~PLfQt=S$3*!qGflnzb z%9y3Fgv6y0KvKq7=3|9JLu@LQQ$%R=kclt1|EAfa{!g0SxSMWU{ngf)Di7OR$iHoH z8J4-+zp?DzV;fdwc(~Qy^BCi9R3+oQ)FvtaSlnEGuGq@R@3HA7PO+Palyfd~s|1~i zrjLkDAi3tal#p?0_wfLZl3j&Z9r(BJvRWtJiO^v1oboY+89k8+evQrieZQ>o^5-kG z@{}m3y>y}Nyr=SdtJ00%;jgkk8>G zX3KhJoXOp$pR6cU*>b>b1u1( z83dX=Y~RbieI7Vt{5Rm{p^8YZ(8&^{1Si->g$!K>06*Fb^F_NVzHP_oW~i{j4uRLCrNI)m&iJJR*fq$OXH1 zw#D}94{8$%X~E43raZMQ-;J{vi_7%$$QUMC$vRoAN6p^mY}C_CaIog;*CLE|GE*?Q zIg2^5;}j=+w&UN+G~SJ!=Go*vlxtR=v}xwSQ*_L~SfmYhlVrd<;$&Ld5D<1Q88#j> zw010uHQ{>d=<8BZqQXv>>a2RP+{80EZT9v0gi4w^gLX-F1HLEt`f#a0w-4QqW^Rb6 zYH5jo&ZXzrdXKgGu(T<|Lx_AT;`2qrNZHpvtTys?P2})};ZwidFo1g?jGRG0iEZ_8 zDX(_p9k~N?&AP0)LE;8S+LP@v`0(U(X51oG`^KluospcOt=eJfq+++eZ9Z%9P8vVR zC~>BA`l{LSG7bEosM?Mj-sXtabPEcEuulyq?bmlH9bQR3rC&Bi5gV&+dIp6eW~#N= zp!X`2!{^lltz=C6S#0hrLdB+)?N3et+7*4d)Y@9g7p}sH zJXg3(AE2|bk(Xl?q#>)NS~a1%zqM|}TAZD}!+o>V!;@M#ml~rQQmhmzn za0pb$qR(7Dh&-gBD>Tzmb1wH2-_VgCTbb;TbNHgxael830%yq$n~E@@YWm;y_5SZ{x`=V>YEpWvscRX9{b^} z^6EKdjG!)XDy?1}3y0^p0y^R-LzPq=BTf;|{cax3E>v9+isRDNmW<35J+w)vLF#mV zijqO-TXwYI%FPrS!e*wQcOEZaAPkr8;)E!Wh*&`$btOzyR8qDOMbr4@732zSzz|0v zW3r`<_A$!|^Mqo69hC9xwuC`B`5D}8O@r+CKo4M2=__JTS+JsrjTd#L#mVV3D3>N_ z*F1QBI@y&vxtLy%0Kb|}5-7Yk)9?sso03n)*bm2W%ptu*X4V} z99EaIT@0S_q|mOQRS%T1N7b;G$8F-{$pn7>>^rRpS>WO9FQXv@+(_JuW##w`Q=UQf z3}cDw$J@(~x9j=SkF&hK?HHlUX3Aq-yR;@Rj1@tk+iV*ymf`UJIdi5wzVux!sEmO3 zGE4K?T4I-zv1;ASg?8w< zZRsNks|teNQ_(D4RechDlb~BT_asGu)ef$v=8C2Z?@l*Qhyy0Y-c*4Ka`(1`6au?E zkD3+`&hd5SuyO%t^P0+(r&tR7$S6E^kMlwpW4D&O-^G*phl zaza!KVHD@=B#_R`+Kx}L4D6H0UnTM4^4!_-m|w$V*Y#Fqmr%xp`uy6xJ~V%M1!dgJ z`o<$8_HpZG0)35rGu5=w46F5WR|jSAPaT{=B{3-}()oUk)3*y+iNK@^Z+M`dM`j9Y z;Le5eUcf09(^_huJX38iRK464|8yN;9{w^(3+UPO7Nk6=3CMt8oc-NA3e5PzIg%t! zG5%?x9n%HQ9HNv|#BHjeDP%_>V&+A!t7DjQZ^cQSz>h4%t<251^*mwpFTEyz3TO^F zjzQus0L|)u_(F{zP|@s&g<_k(Z_*5yG(m%=GZmovG6dkkiece|L_I$-fnz6(OBsw^ zEL*GHI13tE)YG335;L9HE02$h)YU6rx!C&@zT`wreo5eJt?nxc2&e8bV`pzlJBnWx z`J$WP*g}$)ype2W=#zn%3^d#6+9tZ`Fv1qCM1pl)DCz{<-F39Z@u?i?Z;Y0AQfh_u z?P{)QODeFGi-}z`J-zY_Im^RuL4=@|&UzU#X=|$QmyU!KRTn5v*K;-AK(`6jGDke| z5+3vyE-y4UuAtUMnWc9Agl^N4!Ww?J#E3NIVT(M#&ztBFbL0QQ-J+>>PlwPW8nWITWQP^ys8K-F<*q4Cv0)nTBewwsk|p3s{x84cpPB)NiQIV_n` z&oz0w_vqUY1GYm_1B6SowWNl<3RGuB>oBi1AA(To3Iiq2Kfu)K zi36CDm;M!|h^ORPuYvO9aR2Mt`W2BnOGRDcc6{r^!%TNw6wG@)$|kLfF`GL6fwE|* z3F_lP)ghFg;@PcP&`MIXR{>?vT@8Sdw@pv3@P7%-K>Ie5f(ePXP%8iLTK>6RI>ae> z+H|$=2Y-yGkJB#@R#nPSwM?L?MVL5q7Gz~d*vxCIc;qTa5j7<1?1j#-VHzn5y{Jb= zMcw8P%c={;xKA`C=dD`=Tr9K9k?C!|LX69-ne1VGG^W-oqI|C==Val!v8kJjHoI{z z$o3W@58MM5bkFR%`F51<>UH$Jxh;;2!cHYK2P`n3<;Qcu__U6GP272Y8xK|FS(@49 zyFqRth|z(nOn&?p)Z18qo{zmjnK=?ii?5|8=y#Q{JP*9 zXmy4RK<7=jc^k7?EY;MnE$2v8=b?D4J57GsU>q&v+F3u7>abg!@{a+ch*r;_{95dw z=vL;U{lufC{jxwl##!+7FbO`R!MNmP6^`WcebSWqQjJTtWC7$ABpVt}8=fgv5*59Dv12qN=&lG42z*FlFrQg=$jnOFQuTttD>t=^c!k zWt>)TfU!(+jD>KWrB@1TC?Qds>C!=g*u?*3Sd7m4k}It~s9{9AH9_@Hhr3qpx51ST z7<5Xpre%-!|52vV~BLk>!B^Gao zl;%^MqwO_|$D1N>1-Jc!r;!2OA!`##iCx>n@m(^Rc}z8K%J5F3fbQdJ5=4Pe*{|F@f6W0fBk||Y)1dU1SLb553Lly9((yVtm#!fry6GD(&V{rH8HERz> zrZ{#0wMF3Dbp!5lAj={?$Qqt2oZY1_B;N~fW*kSBrXB@z5imn(Q1%P|>1+Lz5x=9a z(TV^xq}U2CznuzH90yBuq^GFn0kS<`K{T`u>fzM;1){E78du3*_1ATCY^|0O%ms&3 zwAc|OGUzXdj|stb90BJzRcbTs{TYx)46Ad)!*>l^HB_Y?po=8>L05=?f(*ahD9*vD zZPl^!AD(LuY+76t24_?vo}y>4UjB=L=qF~?@%bOYcauL=J!u|^ zi}9I=b+w7SFb?iP4gQkP7#e=tqfz;bzwEJ|lY)pWtcZcqNSsix?y9H1?q)$FnsP~G zQ1gk10jJzp6d-#Vjfffj^@IT<0clxnomRwO%?jEz2sWF)1Rl}2(=AM!g892^)ipk+ z;N6xn&*;bCL|mYJMBkRMsEo7#HpdEd3WfVQ*VnT`&SG#^&|;>YKSqRrx7VPma|%2i zjts%kUU)L@9hEe3?81e<$-jmw3ZL%Etg<@1{Y#;2iN5;wUeCyIR=0=$K>XIY4qDxP zV1Vnd+3{7$S@hI;RimG{i;n#EyX5#lklHQ3fYcuC2inbZb}4rN!z88zGQKHKpb{X= zUyM$~CY*!^cU3?3_o6RcBgH>;)-?f9Y(19BhAd_94T5vSQYoN7cj&d7&blT3EjeEv z$I`fEvgEftzER6yI50M8miV$))oW@;J0f_k7X*5P3_J?&Py9Z)cs8A0;0fpwVSJ4o zsqYI~U4q^(6@{Jcf2j1n$dK9Ku-M{V@lK;ytGNEl)8&rMR3gA>wT^BU+-c>>y#vj- zSenJ?I;Sp9HC1j%*Z9>{M~hzLDoxG+pxME<4rnwJC)s6Qi}nokxKDnxPRoKFQSC|0 zb<59TFG-dlpD_JtKFOQt(tjU21>NhJ{1{(Jpb}^3K~(Kn=k;8(yVk&x%0Q+uE|rZy zae0V?VkVDBHNG&C|9C2MGwn#zD~NY&l3mivnTho7{xMK0QbjE4eE$TP-7OQS>O{tm zk?&>4tJq|Np{bvEswIC~yFSF$(FHyg!u*=mnbD^x&H_jT4y zU*SQHCebspEFe$!P8?lwTF|Ll?sD&ZjJgHQ=*HsK1$)R!HI}asu}#U$gwPo~kXv;u zY=qI{$NY`s$YAb6=I4i}n}?Pr$6p!bh7h+^TN|XBI?&C_`Tb-ndn$K)B3AJ)95{pl zZ2baZ6PqM-Rb$=*M873h(6(IEEEdJ&8yC2s8$z;s;#{_tVdOGhVNoT@8w;F%k+?-~ zwVCu;>=q{xr9m~`KrSI9rVB{G8DYSDIseXlXDmHtCE~@-ntnEo6MQ1#VcdvJyzMBM zT=)rDYB-Q1<+t|qWcX8u`sygsV`BRit$Cj2>B?Dm*9Bcm7Ju$f&c~-wD;ZVO%ja5c zH5^;}ACC=cI#Kmb&#$tJH36?^G4-6bvU^(o(#5Ip=f}K~eW7Z0+eLGixb0HtV zZd-KRabR_(hUHc>hJC+Xuh2VGwB0t(oT61Jx9)%764GZS;f$7qQ}lS5888>0xvgoV z)ym_bBJWVD1Yt*fZVpJPfgo)iUsgVOB!*QW%V?MAG26Q+=f^w{_B7PJh!SuYr3S_H z0Sy<|q%x5}$+Iqxp@?{>@mb6?yJpz6Lx4O`8aN-Qv03}5+*j9QQOF|#t~p~?r;!DeI!f#GpFn+#T%D@dF>b~4tyLT$;6>||kTUTBxJDS!hq_z`( z;5e%Fz^j@6f>-NCQFEr&<;twS_~*-Rva;76QhD0_y!@o+oLW9X`S^sgbC5z_vs@MZ zZbY5RkG+gew%1zyYgg5$&5}2fI*1W7P(&n($x;K5fQ3vY5HIdW-=cf*cl+*>>EXM` zu!nR_$}WfErVL)N1c(7ANBtWoFZcn?HRwf@Hrw|gQ0#SQ9zo_&$f!1l5>uY;l4{K{F^g#axtGY`x#_U2eoH&^ck_i(nkrT!@{53Uneu?Du zR67W{rKG}#^}PQI4GK=fr|gJ2mwU~ci~zSg?O<7T5LAGfCLnzN7&<>J1{Bz3#^HVo z{;n-KC}&+y#>%|YK3cm}uiUsgcxE5q3pcnM4K<%Q=eb{mO}X8hVH5QeFDc)fODATf zrrdm}2=%R(=(|9FU~9@F=%s&MWrAEXA@G8fj{?H$ww2OCBW1^`iF<4p>XmiT2{ns0 z{AJ!Ghc#DhbZxBO$lCvlt9eOn={_q#aemW^38FyDZGh}Yf8JRa|2@bKy^}d*G2dsO zb!FQbN~c*)8@oJYU^&F852sY zlIoR-#MxIAkWVtn#|#+k4gJ)+H8tP7ON3epW(rHdC|-sdn#>bt|GSXCZIR2erpB^M zJ1D3uQ!wdeXGc=rYo_9g!4`;DDX$W5R@=D*>u!kL8|)XPi$}_fAr=E=Ds>QnU2lyx z%-Mjt3v;H|WohVa05tT;EU(dZ8I2hh%g~A6ScXJI2dXad27Xaq7rIWq^c3mtj%ecY zvlQuH8`VoZdEoBUrTk0gD#vv$SN$&5x~>JT@1uT`GPA&3afluXXicW#C_OX^z zIz7ZPPu#40{AjRtkq@5nny%g921wIm^4FH?7Z3y(p)c(TfM-`LP@7fj-~q+qmVvxO zBOtTI3ddp6=PZc7k_~~rpdvihEI_48bCv$)e6J`pqi`A26;58ch5U-0p+}qErOq?w zb$_0dwxG9W?0U3xjh+mfGu-QExgushi_|^u8u&49{Tb-qpH6D41qQn%5HS&I31Szc zWORYTY$!&qXsGPr@li1Qxn?<+6LBtq6SVNapkU-jYaQlFTd1hjdbIM6+RHq@K0B4n zb!F_czTY}XB+Is$Uo3I%P zaT4slC|9Qdkwh{0b1Rv1fPp!CMN1ROH9WZ#gX%C3q)5hgfg$}WPjDuDj;-1W&@&(V zBr}rgC+h(XBO~nx#%b~_W9;(QA-4;v=5A}g$%rT4)v%Ol=Zq)+O3SrI)mxxFJ`&nq zLAA_xP%U$iiF8kijV~1h(E-$h-TOx@`BR4lD!|*%g6KLK8{ivTXeKe;(Z^E>4nXo& zu4p>l?CK+I0D5kKCzgDoS3qx_riz6|Rjml1{)sWAvz~YT|BHLN)jW8Ae}oK;m!z=W z$1@H$H4ZI%7}qZzUs}s!Y8RE%Mre6Wq%uj!cu(yxHgRTBwF*?o$aAh$B;?8ECNHnq zTFZAbWsS(F=hku_%i&(jo{W`z>2xM(1MAupLiqC9dM zxv3*RQ8=l3O_ zHUBO;{A=@b!cqgXRa=i&%%upZW$;ORnDS)$0dng-)iwGcv7uvoMALPF{u*S!*`1yA zedP+j9Ny0&;Fl(uN@;-vgymgf_kxo^Xa)DRDtDecT?b9FTd3Zso7zKtcpSSt7!=0` z0{tm~-M&w&{An6mYOYzSVsTT{YOcRBa^}jWVlY{A^u13k5r@5LiW$HF`V&z38nJ0> zk&{;`u0<$B;Uk;D82v7DEUX#=wa9OP zVMtz|G5(WCyl-J@)hG4>R&WAMM1RO|5+_BD{QRyZZSUDDF`(1X`WN zFfN|jZCFhwUR%Ru3)FLTQ0idu7GSk3I|XzpzoxjY?K}I5=L4*ZboX%8{*!e4g5~Dl z#pL&Jfn${CZG<4{#FA6s`ztvaB$YV^Zvm_@RwxDlb>G?sZQ!TPpvE(`(QTuka8EM_ z)c%05>1XXx_xI3T0Q(_N8VEN=<|>(c-Q)mqe$_rBrS5tKX{oE;Rd< zf~u3#0$uZ-L|($^vaKsLu8zf9e%JK~g~GU03Y$86IlJ|OBDJHg`;GZ`!w)qhc)P;b zWko}^_OhP$SqeP!+~MD4=hL3M`^w<2RRk(-nEg-E^Wh0^jlY{t=Bb*#t?}*m2a1Y$ zDFe=&lL8l9+c%${ZQ>-1HlabgA^rzY1|tK00L4HIaw$k6c8^G2-^AjN&%q75nRj-Y zsu|OcOxHyfyM5_T$7jvg!NqbIVl3rz<4lkiWW-a?mj*pTLaqN0U&aLDc zZcW`8NtKRV=nh(6REZ2D=Sy`RHJ0%M7_YBhht0mDl4KBYqA9ulVj#*z#Cu|t5!PIt zrz)R$a&H8JbN>kB9op1UYx%Ff3QKyt)|Q&b7oXKq^V!~YQX5{o2U zqR1Z8?D+EP>vap#69e1BQfbgJedckNgG5c?y`0Q)(tS~9hZxhI#jIK4Urycna0P)^ zaT^aeyBjP^pGvlpPGJ)b1=p&4FxvmI)k`v2yv<^r(RS9sv7S5Y`dL?seuA!(l0^VF z`p#*2vUS#6RtxbjJ-GmP>Bs5yt0-gTB14P$=u9d`)!;4>fbg+OA=bkYV3NYuXIHYo z$8uiD%rypmf>;D6ppq~MGG3ilrlh*ZErA}GI6WPtJyl__4_*NBG3N@+PeN5zb5iVh(Cjb z>M@q@Ja=~xskiMvP7_rGyFjN zyX@{LGWCy0{!zEX?APoLNLRb$E%G{eNlSmyFV1W8pLx2(PaL^`Z?jqweq4(fjT&v( zn?F|WIt9#oEIDL;o~slHc&(db1AvwbmKvGzB3>Hw=4>##C%Moa4d8+YUA&l5fzJNz z3~bm7p4$@O_Y1w-};Q5&ZmA% zzzn<*_tc|dj;}J9z1QnMT`zYX?24>5<>o6e+bD)+8_j2|BU}$W-A_3PkZdmdnNtCa z4#r&qGuq-PoOj^-%bONcv3>DA3z?#BGHG@7b@6mHmh=8&%!KyF;_osHsbE#+ekc=! zxjqcLu%&EhQuBWCdO8+2ly-PYjZ=h2-I69Q8!^2Bh+yhuv+V3}EKmhbbd}ul*t8mc zVoc$=S=>;;yYg|>L|dlc+)1>u_2PuDW*vjE;oi(^5%~mUuEQX6_2$m~4#>*o#%N~% zp)?(ZW1(z%x4iZKPY;9B-3;8LQUgkkTEZnOuQ_{lrPU64=>Sal=Or~UggAXBqXCe> ziDz1li^RrFDOoXUv)@}0Y1eZwlCj{ble4fMSpxUrH~?2zy-ZpuT%@iO*dDmTKq@9| zF+ScNt{ANb^Rl8GsiM|-f?KzDZTNzG55qnT_k0qm-hAI(h;SsiGF~e*8U#!LkOl~; zNyxhb2fl(HKZ*o1Gn1ZUK$p8m>RCjX=kIIP!TluD2c}%$``c2$Z7x3XTvJ2aOp*2+ z>WX$n!9fkZW#Nk3-b95?#$N3L6DKNAY87PR@J=4dLb=99&WW9PWlZOz^^{;nn$$vH z;++HscP}70ZqM)TZd72*U$&cmkvn{)e}%miEooDkME(!gej?BAiLrCU3l0HqD827f zwIIDW{g&Pp5Kvxdo`huQ2Uqe?vEgDZ7HGc8BhSK|Z+8kGZm1mJm~2fLY?Tr4P;`3} zs$2S|Xwuj%oiYxaCX$zHUOUeWTDN=unRiJbD=+a9H|S9}=X)>-g3^+er~NOC`sWDk zS3|i~AUQtDwhScotO|QuRFwjyD{v8m_y>MuFh!s-o%koG0@d(9Kl{5=2>`T%+`8b; z?vA~trxq7sOTuu@h7lEUiLmawD&WOdLt~9dTcr`nmmO`=h4QCH^-ybWnt{QA@)x|j ze}JW`Oh-EhOqgw)VL2M1pzBUp)?$moGW!A@CAi!trrQmgDt}>VSNBXeFA*i#VlYCM z|H*lX!HPX+m}{xIV^KbBtM)2Rh_U}kZNy*lIIk%yOuK-)3h$9Bw2btF)HE|#mq&IxO~1Hm%p2$p`>QTRj1txlqK0y zqjB?#+HeIFW?b-fq+e>U0a>6ar56Qh6RuG1LO{vwu&W)`o0cqKK-)xw{>uDEf2 zJ%yQ*#b1ju`zN$|OqbXaBayhA!PNSAELe*JR(|`~cml#iIk5{|uLk;U1v(ZUQPyVUrId=T z^kG-4LwssquFTDQx*(Xo1FV)E~a2uBSSU6C^IwkT`eZURzDwq1e#1<()ud%Tbf01aTp_}DN^UtJjOp{ zg*<|V;=m(%RmAT70n9mc?gMS#zt6O1M8Fl1wUM`>ID6^xBh#TGp;n&@4X@J(WWmyK z=M10|hme4-@v`7KP)Z8r&k#J=0rAw#FKP}*MDBge1{^_}2I7)f{HTwnNw!qN6*?wX zxS~}Jwe1iibO*RW`HBW9kn7!hg1j!ftA9hNGmOt9*_hJ}9}K0WT%-Mz6g;RI`bWnO zJ_3|x|K70$srS7Yi2u9W^>SqVqZN-CvE&aKh%xR~t!}OBTX6Wad^RpGv$)%^x`?T0 z8cOJ-UP$OWqL?4DGnkp@WQ$&&&-*i zW7$=PK7_>nCM(JOI0ZUuZ6HelRl>m?Y>2%{YiAAt?Et8Bf^#Z!j2I>>wfkibnPHe0 zlolHqto4+uf)o~*(C5Fy&Rqu8+x-Ou<1l~CM1n<6Ja;JekHrsG;4%LOJWnK~H z5IX~|gNZM!TNGpgLJC}aAaCfdA4f)b*!VJnJ3Vf;1E4O99eNCemL1+Ji)pspafC$u z>V;2YYtzqh_O7@+0Enlu)Pi8Mm8Wv=_6#((aL@&cV^bb*$NihlXVLDJt)qA>Ja0~s z-!rX!mx;J2eSNDvZtMDs+6XB%x1;`gXNGMPBW6F31RhjvIr2(89RDK`0hnT{>9O0R zH)L>%@4QBg;mtjzO?xn}I7k&RPvNq(%W$}(vYW70QP$YxbGf=v|tCy+J@>tA% zhsupQL&R3CcLrJrleR zf;n><5Nh&KvE>z7GS!TLRs}bem^gW~b=8&<;{SMsLl)nYxTxdU|TVDZI2a&_-zcp+84F zRe=q8Zb#etY}ey;UZQ>P-HzUP8azT<)WcWwc>6ZsF+9%Vksq*uZRbCTbC9wH3K?UC z-U+jvm7L9AcaqLjxxBOfm9BqUN1Mj0MRA+~yW-LC=>;pY{g>7gz;?kG*UE5QD$B_2@gl z!73=5+8zQ1;{T~JRO^(}B2+l}a)r%U@~yI&10r!@aSyZZ<;l2S+`J2Pk3h9u^XIq; z#lbjTp%OOU_4U5Paiy8Eh`-Yz<1Hw1=?MzX$xY+&0G5xAFHlIZb9z{X$(#j8Gw%ZI z)A^;c&lwufSWdLv1n|RCj~7cTromdadg%=MIov1)U&{7C!R+DMxsf$YpSN3{1vphIu-2jr85z^Grwc z*RKomNWz7smz1~G&X;YCJneR~8Q2Z9{ruTE6*N`Enh80oVIHq3PKW)g9~hS!N$=Vn z$){*v%4p#5D}ai=2CnlpOd{j9zAWc+n=(iud7&4A7;es0U~U*Lc+Z_^YTmm!AU`US zfZd(W9{3XJhSr`vYfVToT%Kx2t6#d9zb{Ux7bn#}uKys<)2=7XF^0Ry;ty>y3NQ4F zQ@~oT+J+bK%Iayr(%(+swU~frkUcG|dXQt{U%yBwbFMzpwH z4~U>4y})ktu>00-KrDK6^`M4pKy@JXX9v~AAwmk-i-X@c?L<3rj17!M?-5LH)N+#R zF1|EMa7bxLJBUo^q2xd`E&}iT0#}+6E#!xNQ~?ifGk*p`?+JukG?ft9uv>pmZ}u(@ zF>YQ023PJEH+B8_yzJ;McgWj&0?!fX|X;$}~EHM`pYp!(BEMK)`Mu5`xp3EnTK! zk~`Jym&Ox)OE4kHTf7c)0|9*?gkO{T8%)oypq*C;Q4@ZU1kc#bZtDVP*ZKd!_lLBH z)N=n1UtF_s_j5DVA;AijasB_yxr98)(^1vAS4UBWi2tj@82*3bFuEja3s07wgbHS) z-#sM=$6YZSdH6l6=3h?d*$=0a-Hn3`-ucfW6f&Ygv}xR36$%ik_6ok|GtH`$`u5-| zo`Z#5GQhi2I|TK;RAwkiFV4kJQS@FFv|utEl0rZ<$x9{Ajr&7)pRK?|-FHAIzZ*)W zmUryY?b~R`@3lxe4RJQE?BKbKnHO9KIQOqH^F;e?w@7Ds;@VB~!=B`tiDlWEbUiF9 zm8`go9}N>72SN{>d84S*nZ$!4ozfH86*hmr#dLOBdDq2k-}`a+QD?UfZ;OhVu%rh( z6TOd4ua`&fyooM+MS+nEcJTWE6-3-?tRPt`$)Db>s}ee6kDSW|8k#D&7d4mw!&bq% z%Vp169+M;h$da_6cspY4Z^BWTs05I`yHMkq_+lg$*u*7?W0 za{b<=iZN^bBU3M0>H}O*dCK*eFvo8Wk)P`KKeLF5a}a(O@r*l`&v7bg^ z#%0wl71@!4axt>(+hJW^bDO-%av?WL`W(wKUlr|#ABj8K5))AV~49(9s&ngtzr((oejFHV&b~#+O3p9{+A^; z(b(ByD)Z7l{~>(9NY%s&F0eUUWlV-uyU)M7v&{3m7Uoa-a?t^A4K2C*@WFQHMVU{V zOAF&hPuWFa5p^0nafR%lp)W#Y&9MfesFFRchvTk|zpXDxsT&|d*Qh7d1mRg&J57~^ zJ74@7IdJ6xwl1BHq0CKjESRnLFmX$)-2Uw3UD?13fFhxE^S#DtL3?q2+LiEMc-m)! zHuN@z%3%P3>}N1vv!@|-bsGk8CJ$Y+m<+vz0`ETZ%Mp zBkX1mY^It`f1M@)46nvU!(|rgguypN^a_Knqqx?sHea> zjXmKRpZCH&)a>f`a1|f?zaCb(w=}qIxaGx94%mQfzCvG8N+Xt(!ij^CyUzg~yC~$g zj}I#L%RtZ<+h5{?A@uRsa@MciX8pJ!O~u$lC+#ANqqjZYQzN!YY0zH0fHx1J_f)nG z_zQNv{RSZh{m85mE~)N`f6$|~9ffDQSO?1Y)a5kRY;YMdg)&Qkq~?_8SL}7i`umfu zX9cp6tG-5;n?)G%7Y3srt4doZm8iUDsCZ1A%I~I~%I|LkxEvlEf`Jwstgvi+^LAuP zoU%0a6Ed34v@ZTe4E9Ca`g*sJ{0MvN62|nyySWS4(uOp;P232S*$JklxIt54L{^hq zf@*Ia%Dg@u?@epqILYE7gzQzAyO$?umn}ko!!m*FrvFuFgb6JV~jt4lqfP>m2=AtB1Ez8Vo+`WE` z&{#B^c=js$nI{Psfq=Qnp2>X)oUjA%^FAzyMT1@qYy>uXMLzljo}5^>zPJ;jxsaGJ zVsyB}fSd z`yg*AuMI8Nds>yqu1v!h-x5Jz0cS3XTx?zQw7n0ZSt_1-&^S{Dwz??FdGs3h-`E9T zPygL#LbRQ_KOJENkWnTa;t_4{-n(=BZ;sVEQ%)IwAsNcCyW-!o`u9XHaBu z0qfT*tG8;Gd4O&-f{<2UCQLxi#y#5&vywX^obHPl1naTyu_zbP;N6#OuQD6Y8#axM zTo;26)tmbCKm8kF@4E7b!&>{mu^zYh`U`%uoHF7OQyFx0t2&*u)Eh!gOs50$t$5!+`K@qF97|f&dVc0 zJ4?>iR^*y?I?OSl$>5hN$?-%dRwE^3H2sVSemeF*$2IR=v5$?(^mQBSX6=bC-k`@VuibfLA6E&lvo6}F+IJT;~mz`#2rxK{Y^kH$F@$@Ht9j?{F80{iNx9O;=toe zd?wux_A`p7lJdlgB@-0Gq@}lK>+j`he`QnuR71qt1d)s?cy&nYnY3wcc4o(%iGS$M z^Y)*6-OD|GK#m+!Kd4_8RsC2NmuA5G03;#y4vN)jWx0uKIBurJZaBiNC?Dbzzg$5;9rKep$*_3CF8u?-5lwJi9-GEt%in;{qho zN$d^?6nX06pbz%Gz^R3x@w!u(3$OZfssx?W-nM zS0wj9b+p~6!O2j2FsJh<^AcHwrAL{$!9RPHHR()ekXCA)=Na9+D&sd++MCH(5^Fr! zpKgI+B?r)lc$=9GcP8fIS6asA8Nz|{g7b_z%A%{$+sXTs*`{!dDc(MEr)1y)nU}=4 zJwHBIOm#1|qK+0&=hc&0!KAlF0ZZI1Pi&H87LQJp*qqE^K1|EwlCDyD9K|(d5q5hg z{fg-Fq}x@Oes5@(Ck*DK*od`Y2IQJJQT2cfkS_(Tio5G(@`u+IzQiqO;>bT_il@BS zsR@9o4zNs#H1-tD2~_`pj6xpeoZvg%CSt;DV@>5LGS*XQdB;-srW!3jE6JFhU#T&| zj3zHPNJ_872OfS%p;H4_Vk{O@Rb|s<7UH-zjZOqsJK*es3=U*))egYWne=#}YcJ)n zkdcXoy#oO0Z;m$@sAN%0eI`)uFN$NNn|2dFm?{6oj-{ z8=fZ$c6q?agx7flX8xcLmUZvh8|c(RgG})~fI7{=PV*bNt0^9bAbTJss?>F+f?~N?vh zeZTx#E%yB7Djo<8x(e*jVC-(Fn^`6wopkfL^O~fhwKMHaqczYoWyi9=aaR845AE27 zktZ*9`dJ=)2w+kofS1pCphNvpS$yAQ;UC3LoJ`R^IZ%4`;G9*CLW*|D-F3X;#Hx)v zW{{N>+u2#kx=~jlhf#;|`QZc!2w^-8Mn~?5@-}GILL(!0d2R~HBOUHKs9?~W8#L(( zjM-afE0#4?yMzZQm~H6>sF&e?I-*v}YkT%OyDOPhe93h$5KVJDA`n<0C`bao<>UE4 z?zg%=@0%Hq_jK{P_0~}9nh6tUoc?rASQ51Lnno4AXzsnd=rl4r2#rvrs~Sg95pj;J zy#u6+zq6O0?Q7W0ZEYf6alvvOFdZ}$67DD{Wsd_oPNt(NvNjUr`OFL>TR_Mg(B-bJ zt8_i>=Hsg+(SSf!P`LPbiolr!7oNSwE?eW;drqP)Ig;G3Dk!x-eJcq$><}54ie|Ew z<2M7&$clno0!RIHXN|3tJc#I9!))47^jno#|BFS!pfXsUm(7yUi>XmJ9gLnA1$<*Q z{8JC1J-<2FOQ!mlE$n_fPz`(E2{rC+60~H)Sw4*x_l|O#Pde?=|f7`6pE%=$2}w$*|G4NSp({wau!#o z?kGlMFf#8C&S`aMr+3Y1iEEukv^arS!5$SCXs;6A#E*XLCg%$%Cgz7aowLq0xoMT7 z0X2eQBfCcEcyW^Z0kp36j?4;N8gVW6D{yC=hgLfTA5^`i==&@x--Ww1Uer9qwdec7 zOMJ|T$HX9W++Xj@H%FoiW%Zohta@*y64HvTL(Or26%3&H{J=W(9~~LU`g>j#J1vY6 zKYmVV2yCD;FVKUxEUe%1?hB7w_F)dzEd3-2d`N!v-U(HJE``BA@eqF;lUK75GDuiW zHo6Ni^xZYUW(-CTf{YJpi`Af$mMKE|r5Gww-_4?6!Okj$wmM0S@L>OsLhY#ER8$OdP(`yuT{*w$!wpK#Dtu`261#C&9>0*q^(xFHh0#?)?F zcxll^@QbFkh1S{Om4K1(JoIjz$M`?_kcSer69o&q3jTP10$Dg0*~_DvL~Xj+8gvZp zb~4kO_w1jAq%5Dj4;cy(nBj+eP3n)FSAs|G_mglRpS}^Zfj`qP+I2f^N$wpq5O^{< z>pod9eKHA@xUYQS)(FU4EHH6D9HtKst1aqYa_!F!RP<@GR3U1+qfku#ap?2K7Nubj zXd5X;OK$j@k$=-U5$@*E9YdXvy#Tyv-&~Fm5ng z4YW-ynAR8BRS%D_@%TO42Bl_x6@|FI>%O$dH%!h~9XxwIwEm5kKEfom!RQ((9*O6Q{6UJ? zxV7$nx^Iq4wVf^reywP!+jk)aMBi!@qmO&sl(!sTk+*-=CrV#;X8eo-bb2i6JZ-xy z-{#Elyj)8=tV6MTe_D}SI$S|Ts5f;mD|{pC;_SsR4*Wo1wg6!Tw_JkyJI0GG1R)@G zlNEOG3$r0(ZQ`qwJmPE^AVrai=n!qkDbEpRj|#&j z@GFM&N>!5D!pehkGHcW_P~$+{qt902Rp;N5OIv1}1L@D_8Ir684m#rr zYetus$}nU#Sy1|ME4RY~hAf$bzWz=4vcZ8^R(!qB_^ioS`Zq_#ZeOH?PZ9?zR^<>aiOE2Lgol%OXoZSoG2eHr-l=9I5=&6gv>femFjSpiu3ov&CkP+#_Xy*B^G{kjIlUVc86LT9i1lxft|C}uj$&cappN%gYu0!83 z+CxWRlf0HI?Gya1y({zy;eeR;ZyZg%Aa2AyO^+7ps<%|9Dbgco*u}5UELcxDbi!>s zWCKv1i)}wgZiKdi=$A1&eFJZ%66}VEw8r+_Z?lW7w-a?}ynJ+y(4H5D*`5<9&gC=PU-FWSOWb_1^%r2RlTMdEQ1TOZ+>MBSSg4QYY1bCy3$ zag{OaNoo3&5nRAvRVZ%dW+-d+ie;YO;0pyYmKjBaE$evJ z2?>GN7au_SoXV?=6oX{N3vxW|l_S8#=@mVl6=o;t8BYK1wFLQb?dVKie?_tGF# zU}jdm6~6muB?_5M8M)Tk2k^J@Ugv30nCnWHS-;>)9c-M9X*;3S9>hL!ZlQ$PKqd#S z@bdz>n1JVjjGBi8ZSdosO7`pRtivQKL?7Nzzu8Kq&sXPMy34cXNafSJW*QRv#XqQ9 z{ptJT{l`eu2EoYNC&P9hSrKK9l3w0%e<#4M0tM?Tdqpy~Dw=dny}1<7+lATiCC8VW zH-pqx7!U&rUJXbJs_{klbREeh4qb@}zuYU^b+H!XBUS-4rFX`kAtiQKHAV>uj3Z~o zw=QMnHLgWd{0qCL>q9zZBx$@t0dD(eVqHHoWxE0W%VO8g@~SROxbVJ6;t-UYi;qLn@+6LY{7O4#q4j`op>P`9=jns!V6@;b26ip#?{TeHS(7jdJt)$fm|fuOIJ zN1oqLx%%e>*ZtwLo@v~@lzxoBaf7B*lNb1N`9K!R?3pjkSBaKDy%{sxaq|GTS>?cH zz(dNSpomuU&b*l`K>=`rzh#tC%*x8dawCyNj7*!r^gz*aIw0b8jsU@-1Eu8ZCJWkL#D9X8oz-y2JL#Lr}uw-U{5(A8NipE^sx9*vI?S=qm@cp!y2xI=OQ2%DuF9 z6|q~l=XxD<@Bun!0v-st%72B2U)e#&!O)+vhGtm72I}gy8bXCU1+j=V!kPzHJ7J)aRx0Pg}GUx$&OMxXyxV& zyXwNkC;F(TX+BzTD!=~gQ|Z39G3~^^zheVOfIe>K*day6u+vykbDSGxA3a7_)45z# z(Cqog<0yfTH3ZSb7ECU_h#6hPJS^zN%c2-Dr$d!iUU5Cc9D_+s*2gg9Ht@o#;c4&G zaT?xGLS7rVgP^^N5RV_;VasRA-;h7-fUfO3&akwc3G8Z+WK~WqSIFX5H)$>As~qie zk+fhSGekCCmCcrY^%)or)Ld^iGUSjh%ayw-cQ2g0tHeZQ18;!wa`(r>7u^5bi1VUh zH<2rIxd)d}KX8OGem>}QX?icN7!aG9*5ZXhRl>(UWQwz_Z|%;6^@9T@_2BTicX;lund4*dBJ{nvyd9P}?) z5FPEJ?>f|}@>Y3iy?|+`r%n;o&9_cG^*)pgWGRlqM!o-O)6bC ze99PLe?g}?4CLS#cBhhzbP}~&Ja-rT=G_du&-Q#t(Ohg)d;N;4`DgFtz5RVefi1Zk zMAZH@O%e%q#*5;t^J|x`sV@$+x8UVtyd@adH1_q{&ge^J$a0(Q+7D`Bt0Bek0V@bJT)#g!ZlKWEzW2d&sc9X? z$-mp9@%sUG7W909$vu|RyvetYx3&BAx!&&8bg<=<=|G?M*iQWDR%;lSdH$|qsl#0# zS1K;-D%4mp2>cns$%y>?B(F%JT+V^sP1{k81I3QD)pV=uCG*oWSaPbSdBf#XUuo*$1SipPe&fyTgz+N z%~z4oR^1o*5{0+c6ek6M3?G@5P0c{P{MXf%ssXm8*krfs@gpO_|Y7e(Mh$ikKF-yl? z#IIhu_WAk^Ao;23SC(jcZYnl>XPCAtO}M*p`rAj7l>^Oxm?ezR_pCzC5jpDotW;fF zsG;q};x?2+bio}goh`%~)XQcE0=#D*qb1buHIDf@gKPe)X<}E6xX3OV3+FQ^`Y5Vd zaaaTK7oyD#gMH18NX*MpXoGUJHRwEokPk{|p_4pSx*o5>!@W)|peiVRw( zYrf<9wx#goIngtv{Zdr8ColyAtw_|_U}>f8NBnzz1Lk)yp9HS7Sxf}Xs&lQ-aWBW1 zqKjEjUw9_J(Oc=H*pqM!L5eepy3$!TI*=7O#Px12or{7??F{Tp?i$mfEH25Ey&fyW zkz^REWKWWe_e|0F9E0@e#b?K`aXVfvy`7+r)t9~LrS}N=GP%hU{rK^@E24}Rav!Tw z)k8ydp`2@?YUOs>)V&F^*Yh!ogG_a{DEXks#^>WnVl9hV6M+mv6?e*^JQeD5Ii>6d zNWQ+|qu>HFFDq(4zrKGVVwy-J6!453T%G+d#p#u(W+KY<7)i(~jOn5>+wVmEsFpX` zu9Jq0`V_04oql{iSgjk$2l_RM{6!Z%oNaqDN65|lN)TnTw?|{lC^S9;r{dNfm_q?IQYqy_i*=9LstVCRvRb3`9`%>^zu zk&)36Q$3I_KXHPzYzS-I&PB$NnwOfmT%L|R=9kOl!sj$e$*CN+IB-m=4lrF)j!a7S0knmuPZyw|>#%D~dc_`NO4sXvdpZ?A`lEcmuyZTt7|B z5;)V~lcP+~HtL1OTK^cfVH-%b3fY#sHTy1|I4W6k*RMuIb(_`9mlV$dMcYo_yE4l9 z0o*mwAvFL{8;m~pP`=&D4{n|w&=}U~28|`2d+GR!X{TI=j@Gkok0Q)*fEnK3lZj|Y z0+8XofhU5n*>1SgsZ(u#4ov=cuwle|Qd>Pe`!M@)LfyNq9+n1E%o7&Re#1lj85@9> z6~yHW#BpFis#gd_%%JdQ0pWt(cfw+zD=2{?iTh%u!iR&#x-KHP@03PF6choD zt@@{)O{g!M8qH@mT=%9!W?On;g2mVnZZU`i%0G_E&c_Oui?!8=wsRw{*tvA z;Uz!$y#l_m2P&qg{m9%cF3LFz0%xr1Faje$D9y0a`>82+8T()zZ=MAT*MVW&Qk~Y~ zO4MepP*{6Od}=S0p<-Xdq+A+Z>Dp|er5mt3X-t32pSK2u7t0#ykQgGECNb%GY289C zNA|v0N9*C5ySBFuH{c~d#Ruj)e%qQ(fpbAKTsls`QtdMXC34UY(KoixWp_-D0;qh4 z>uE&C5-?h>QEaJ($}>$0RVT0a%>DFTCTM-K!G`ROc9&pxp6;(Ly7by8*!J*)yEH_R z+BuRQWI)&;fwre9ZP~?rg*^K72b#n9_?@Wcw^r2tlPkEn!%w09A8rP&o0$;Ems3fHA zL(4uMN@%BU_813TG=TCDm1)4uy!R>;0Kkd?<{CTXz@{`pc#?19!@dK*Fw>Cf#kXgp za;R=yaxA)`wyK>|W3?|TV;ye-9Nk04io5{iH zNc-YlV_LA~{L|x%w>Yo)9k$nr0RD>jnUyzr*lBtY{Me#msC zHJpjU1@u&38uhgQ()y29^Y=05fB1SGMC~{c^0z1Z4jw+R^)FscAsS=3*grwGkzf+2 zp6$n5dCz%nnGNLm#F;#e#Un$UU`7S^kkWZ&f>)Rv#Em znA{6XoWXmS>S)Mi`{QBz_lRG*r(yR}Sz7q**>qCT=Y=hlg&lff&`1|v!%V)`sUIz> zsvJjlt44iZpL%?(SmDM|wejv-FyzSgP-){$lAweP);_fbJoS|JK_h~btO|IN9f1v( z?gFmSARVSba2@2UJBwZyoTJ?;ViA&(?vBORdFxXTyEt)FV1JibX51=Ke{z zK0!Nm%0tBU{4GBp~n9l~5CzxS|+>9YcpQPjvakbv7Ee*h7v4Ri7)L zE)3G*v$^9&kZN>chiSOlJ z(b)Mhp_an6Lv_aAY}crzUqA0S^?L@cj23cR{p05UR4Macdg3`~KT-8B$k7Y+rGgr& zeN%G3VkaJl`!p;4XTN5hjpz`~R$|}OyYQ86(qkuy7QTVuZ?RFk-mo~K9pUl&)B5xW z<0F5OQ$hR`#068cOB3Jw!Tgh20JQK@=CB5fDKss-$hd*DV&J95)&>c&I5+!Anf$q@ z4%xrP*E0)?RhIvnLwyz^sLog34p_njqbuR~5UJR|UDRj3O4WCqTP>t|ixPAeg}!Cs ztRAZtUisaY-xF8av-rKGSIV+@#9nSZI4f4~?LT9GzkC6u*`hQ`-F0B#h0F^RC&qE zCF@6YHy_gjBKL#xwHc6nH^((zOSW8IJBq5u6*`@h8>^beHl*Ra7lK{BaShwO^u;@y z-2|11BUH*F{8tI)UP?c}#{W2DO63T3G<`ZU0SaL+FuP-<{xv;ddW(;4Duv`u@uLd=ER$dYXUr%bw!6FM z#l6wBDt*|X)pp-aisFXw%8^bsFWin`kja zFoW1$VBpc6*7f0b-2jPZ1o~l5e-n-i==@2Ap<|h+v?JvQDBGJ7wDq9_SVwDNGsD0e z1CjCpc=>*G>U3x*KYCOu*+UgoqJ?fyE^7J#>phwtuZZJc#OK%(!e#OWpxBWAwa`87 zal7u?&yZOm{lGrfz($q*SQ7o|xTg+uLf~umB-;XutbgI2U%9^>-zqhQr z@NGo@qpmrYTy04tR7%y{4jQ5eQUrO})e#Xw?&PFq7i@MeYvp_% zG9_LuOI};=tS?D)B7CmDa@-`P>v-o`sjr+zk1X3T0X?^#&e#z+ zHkz1Q9SyPvg$QJCcpz}eA(`B%xA?>Rnq zUzOPxYDAe}`nH^9&^%T27x~m%*v1hp#x@K#Uqr6F?Df3z)}^FzyMdfU1c^9S=k;2y zV`RPRLccj4f$9yqy2y*-$BSw6ycTnZjKuq*-1c+L`}1pa_ z{#AzZUCXXOxKGsR4M{c!KHuyNCqnNac=_me(;Evci8cq>)c zgM@v>R7^o+rs3lO)D7CAY#>WU8G6)mX~@tYmuNj)oQ_BbcbNtQ3<9%Va8w{m!&=He z*fJ9V@?5z7uTQfIL~E0wsX}LSFuuv*-) z!U)up^1nf^LZS=S=cEE=CIxvC>l&Im>T&LvDukheoPh?AUVW3uoSLu_`LUqNn< z57H+lLoXszzs==;{I`?qcpEw|ZE5uxU+iF@gKKasWc;EH%L z!FA+p^vehQUX*YNkncupc@ecE2@sKBchUz0%B-D*#i4_%?~>+x4!L_5%nAqO0P8Ly zLAyjwU*z4=?3#2TLwB{SV{YYs-MQXbtLVwUvkm{&I_!vzXX7N$pW21@Lca_>hAaj| zi7o&B7&eGdW@x;HkK;V$IoQ;G*cJEuh8A%>*8jdG*^wp5l(FW|cO~MoGR+3jIA%=f zA~Dd7)gb?ylkq3d5p@c+1u``|=!}W2E z&SVZKB6r_Lb=9Y+m_F&n5tuc+A}W55YH=Uk#dsuT-ty8x&!+NkJ!uP@Au8M3)FNJD zgq{o#)JGJ;L3|T&57v)ikEZ3LpDOrOb{xjQ=n6?QuaF+26Ow5D;ERULO23s#9x9@W z6-fvJkZ8{Ksp=7}XVrfe#(zCO2+I=GFa9KygePieASU+NG?7NxnHP_X(e+>bgOBKp zHGrAp);$!ai@u^*VNP~3-CuxCbY40Up`Us4pYF|ATfjd2)f4|!@YAo8Q&3WW^p2+? zUW(Z1)7g-gL{;g{s3_B_yl}R~Qjv!)%vY$N9Y{nemK>V0>*-d#0)0O*jGpU#t|!x` zZPNu0aXVmNswb%2p57F)eJEM~q?>n{Moj_;o2B6WTORQg z>K+Q)GIku&tVY?!J(`c$IQAWYIoyYgS4CNK*m@Z|Tyuezbvkcg;SlI%)SXeIZ!)dG zByLR(-X*W*qA9v`wRuB*<9@#Y(%rkdzIL5k!o$O`;Rl~@qQ4Lr8g4>Kf_nU|6Occ| zZ-igq49ts0@=XiE3}S~k=gKD^BCTg5rf!AaT_&0ereiXTLKTYJmrrmmQT*eBs(=N{ zS#S7!rzw^Ld1*@S%AYo0{1|f_}MkE`a1qV|FyXK#v<>3Fr z@jV`zAD}$B+C)P?TEiJCB5_bU@iCfW*>j(2>GYYIjPX2I zp1OI{1~(qC@N&o}%b%CZyEZDd6;woT+j~nnEt+>`}&(DvQe>N3}S-Jd`~QLwNA?;zn$nh)*pL`!pd{ zQ0o8()U&o5Xwciz>-hd8NkoIXwmp^P*H{!ZguUpLq6h#IT~Hc&T(uBXuvUHR;xj6v zHvINUb(c!>rwid*D!*CTwJqafrfHi zI?qJC_~|>m7Bj_y-p`-eG#>9_#r;sj74QaF&o8Nrz0ze1)Rj1LA7g~#PQooeo9F>z z@ocq&tz)*@s6od$Yz&T{_M zjezk9FR_S5nWYq5iI~WYU~EcJK%Xc%3eQ{+qznK{`xmam=#91svMY_r3m1*5_p^&- z)*=VxEPCC+j2nE^%k!TN0-JZLoNM_R_s-$OpYU@zwhG3gC8~ZeES_lr9MbSPW89ry z)I{I9Hnb>NDsz)$eVx$G=S>}YG}VgPIE2tH!SQ}4_I^fcban*!^1&KaqP@5R41sr< z*sRWUr4hqg#Sk`M>0kh&IV0c6rV6(x0BUWN2mL&It@Nbnx$h%sm+>Qxw?gH**G9UY zM%x$NttzT3WW{&V$ppJO8IEHalz!*Fl1|>`98dGox2irP7M3h(0MkO1fe#qBU`;MI zKT2T@C6|e9ry4LO0*28t0n>qgLl%C5AviI@rP;@}MZOFrHW(h-MSk({&}>4jRXW)O4vzop+Jz3)@~{`}#XS zHQ%~~bAw_A92S)uS{83?c>a9eVsM4j?i(`F(E~tr@CLU0(P^%Tl4i>)n*D;+dZN}c zKqX+=D}JmtqcVVp7J5=gCEA}WDBVOQDK!{#WN8JOy&N#MG7~O@ivTqEIfb7nGcTU$$pA)pLi%o7NLbwIghbA!xAPFAwpIQneI`PA+2e_i| zg@vX_%6(oMwAEM%^o(c8osU?~a~1=|1&XIf62aU;tQWK0WeVpAgx>H=yh_{2#>+eX zSkP@zPe>A0lQOmlH7)Cy%Q`l>$#|Ouv9ncu2wTIIyf&i;9haXZbrODix_g&wAKq6) z=lr45pjy=L?I+yC*%1rij*)fAofcqav~^XX)n?c4-%CovKL_c9uLj|LWe?SI7Xe8r zWuIZqapqOQtDnWX*NTES5ii)#r5Gk^>p!+e_LbXiAe_RD&F=9w?kCQCmoCmgMx^?N zFHufU>`5{!+Pd_RbID6>(NQ!!0S-VJi=8u@BW{Coh=!Eq7~}L#uXU)6AC5q9^B0EW z{5q!1Jon+Z4`^HYCV&c<{zym?G>9gZ%Qla!8~g5nQk%&S0dTYbeRwX|K3ByaPr(=p z6eboPZw9+5#Nv66xW(@d%7@jl0MvZOoe6R^)AB(Pr~v)s=s*p;iiBc^n(sN{I#c3E z7yTb>QT{USYUgml0sZ3csY-`S_vm4JvK%+|Djx7oLNU8w-5lqTDFCJgr?YU2% z)~~hT*t|c=d~!-T%Uw>qwm<=W?<+G(%SVqz@O49siPQd@V2n~*IR3NVBp|t z?ZfDEe!(j6aR#kAf(MaN>{!f>uo_Nf1a!aJ5Jxxm`^e--dRAE;tx z;qXntOvx6QfB|_6EDZCtX{X>rR!zBsoE2Cr(Ab1+MEv~Ymi+v`x+Ux71qa8@xn(~@ z_)s7(hWNdL6EHm6j+$mrEjl?YkUN)w9M@Fn4XvRyyEmomBYh%CY(i`h0qJ#%Fwd!{y` zo0Z-7#NBA}+z?7kFYKQ=_hc)@m=yA%aQIX_oZr?gb)K)#JvuO`IM#!yb`{*<41zPdE4d8DeKR3pA?F`nE znqO`5Ju||UL)h;yy1tdcw7Qb0d0tEDB+sFV5%1j4+sCeCos+O^`OF2)S{?>JTe)46 z%Byh{G2Z4^CJ$62wtAoeaYsuiZ`Wb|_~tmgKjKRxej7fdh)SO67k+m?TFdOqL1O-aAI?JQZ-)~CuD@vh0PYnXoF-t4;HaX&UH z-Mhgh*Op5}xrgtOagi!b~!4f!e`>;Ek(41zk#2&B3xwQ6R4A}8z ziX%nflbJ?Gx#+2N7x=d#oXQd9`-B6>3zdgcP3>mv@5WXkL}Rc5EkkSMSY(X|&U71F zb@RuG5KY(wQ}BRBf6=@uFOV8;aXPZRRnl!4T>QZ1#oysQ+ZLY1p3Tka4AlIaZ z4|sD$w6EPR?5y&f1mfyw-?v+%GX_~5MnkLfLQGL2jWbnkd+K*s1Lm19N2!JDf4*bi zfBkSn$xy^)4k$rPd0mtoETMU7(d(u^HHzGAN9{gPfvM!l<=E5L@Koly(}1_>hz-WM zTX2dM|C-xfy$va8A(Jr(W%G`Mq`2hkP)&iqo-_9~x(BwWK@^gS8unk1mGAW>0} z3>c|&P(x8glu%&gvMZu{>#Zon@`OKy{sx2}5q0qV8xQv(L#8=G`bgt$^iA54I1Inq zT|Pa!Fv6ms=<)$yCqser{ zg!4AWca{y7_m#@|*9I=u&{J||U$Q1$jb_zYB&G4WFIjkMX{QqFN}wl+O7|?A8VF~# z!?i3nkWmhjSItWxa5aBMAsEJSMJfF=L4kwjMv8dLsW{K^stbkTqu;6I@8C=oUH%3C zPQaOCrbT($9-XXux~)e7T~;ueVLDj~VytC;GS#v2b&rlGYS&zylfHN$yss0LD(>qR zGFa(>UCIMzbM_(o&o6V7Kf1Q;$UVMebfZ4`g#M=!XL4D+o>C_Kz`C9`7Q*EUq=v?s zZhCGOnu{u^N<=mz2%>O|nZG?NKqi@J>Y}+XKl^V#qc_^|>GAv48L_eG*5f>t#LVh9 zDd{|wvmGfY9oHPT5GH6TRrtgoPSl)mj?;hsEbqkZfBB()cmq2xs?JlSQ9wH=*1RDtm>xqQi~*I?EEL<5PM?#0Yu}AIr=!#aJ*YSzK+W*Ca z($ioQn|TjBM&(mO5)Rpo9tgQ{{1c}?oK~0z0f!snvT%3KW|Oe z5!OweSE#>Ri`Zai-|AkWxx@X&>jZ!r+v;o(!%DDCrxH?(NgXrcpk_Z;1%{>XMK)98 zPuXS9f|iEEA7&B|CIGK6dl3D2@t9Cm7Fw&a!j>zav*E1&%M%jK7J~inMHa*}RRQ(a z(F3t{{V(n~#FDWi`%^pr<#%XDtwNe@g-@OedzG}=y1aSJkYqupfrC@F_B@szy2=jD ze)83-Um9!jR8M-(Bh4`BnmEbDpP|Rn9Z(DR_hgcq<8xG)tK%B_)N${Xp4`b&6EBhb z3PGvMkrCP#RH8{r_(tTewA^>wD6%{^f_N!iFrCPY8p-OKWEvsztK^5lq`o*aQz;y6;9cgd=%h z;65wf0FBQJI#b-o!dtExA}4jNeq7HNzn)aA|7u{?1wLWL%axxlGt--=wYJlaOehb@ zcxm34_N6J)k3=*8?V7&)@y%6)+sL z%7Wn-v{Ah1W~3BdMj{}R@6r*?Ke1?#G$$x&w0N9v{2y3pjRH?(w<%NC*c4 zw%O0fHTr^NW176QhwaGs-M`7s8HA1f+q*=+%RWr33o36#-lNBoexQ3@>aub^=n)({ z6zs|;0$Q6a7NSg0%*gKgqBnV=__lZix9Vx^lzrz@_rmJEz@&AB{z`M_?R-a`>cRtG zEP16V2fCVg^`Vrx+65{BbWyYWaB^)Yu?E|Z9uv>5eqpuB4#@($HGf?I@SbD?cq3~; zMZbZVf>f3;RRcT^8G^@0I#d}g(W|Lr!48~3vtS;y zcsK^!k-FkEtwixLBP1z~X5zJx+orY$>qJKr)JM>ZJ~Le?KkUmj9SQI0oRqo!J!zzS zCMn@gh-8(7JU>kEy$!dALs4o^KaMqUL69h;aRYm=3Z!z}K}rvf6jZ_5PnI4P_Z(fn zs>9mYRa#LX_V(WVDawsPcz5l8z4I3u9yI|lDsXd53=nYH|a^KjlG03){1i?!|4P+FOZd}5&%x+fA`oscXd~BBCM|jhYR1TdvO{~Be5+1inrr>t0id*t;Yci4H^%82D@WxyEp+bj+9)DbNPa09iI2Y}MYLvGX zc$`;5loN*e&9`i|s6@;|0B>oPbh^xNtHpzcG%_M!>SRjskI1HKnysjHo>5`E%Asn` z-^6E^npz6^lnbaT%9M!zkDnTWPw@Ti4AGID7b=JlBxRf5T)sXjfVh0XM0fp6 z{`9RIlK^KL*I-Y_mU>43$m%*hQ*SuFYcmf{m18LMq*P#e#8Q_CG#kne&iYB(h!VYA zJ{B69z@~H1F!0Ly&E$*`G zc^piqF8gmPZ$-4OBmH?{9u+VwN!+|-@=8kt_3-r+Vx;$RvyeHm)o{MJ!y2(c9tQ(P;3q2M+!P8v0JUpOnM_FIwW<(tLKGn<;O1X zc}AC;u}i+CB8s=SUivB0#}`IZW8-7rb0dAjn}V9N8A_$LrY)LB9+)82TZ^WruEsdw z)JY2x_Bv9CGG^aMSx#Q)l7^~IVeUQtbi3D zt0zglw}Dh(BOZHRzs)IG1lt?j`4gFhfD&KGww-+$Sg?Zj3PU2@4W}>*%XO+s4Hdbi zfKfq-xkx%6*;eIaU;A{S%oc++&tB85f=#03Js%ka<#wF2S?lviFj|8*RPqg_#Z+3uGdV<-1qsgq{ z3&1<Gf;6V8##l+q>z(|OA*)+UiuxPFwS2&8t#mp5~QFeeEH>!BNZ>v|&8b1Y2c z<;o`bqD^>F-xlHtdzWali(Vjf9URuy3%;OcILqL~>(JDH^+GlUabHCgyWQv}?5<&oN+a|herafZB2ZYEaJIY&ku=yaW5;{K;*(J(1vL52Reqm( zk=4#0az6h2%P(X>6c1+5dDpYEGRKc^AB)|?J-8GhN4QwN&2A~dJGc2n zG9Glo{hu8tlQ%rio1cr|U2b6r0nR&Ml0uwS*a5^7eQ$E&PUMCKy+nfc zojq@#BYM=0DA1#J4h`X<39&=Do=5*{4J1d7+PeL^EdUy8fdH)k^P3q_=RpQ~joT0| zBvfhpfpn5G?W3TE6qV5T2RUk;wmT+`?9mzGT){2XbqG3KyGLiH zz>a~7|6URJ*cJAbZC2O{G`T0m*6N6-Gfb)4QHTK!QU|Bd8}{!fL4OCSBmLR^*ne`? z{{Kh)#~(F$3smPpIJ=$I(|i;;zjK-EW)}SZPjOcs5B1*t6`>;9M7CU_EEQ@jHKSCR zB9S$_5yBX<@5)wL<|0ap%94G_nth)TBFtp$l4Y`+$vS?YFYfX@zx({2`};l5^Sk|V zr`K)1-#MT2Im`QV&gYy{;K;Ht=l*@gHKemcFE}&XwYikC0!rpHwddWYRPL*~lND@Y zBNQZ4lTbCX7RKI{t{K(r&|KKZR#Od*l^mec-kpcaa}834Wl-_Gtzd`QgLX$Ng5VcK zTg^L=gk@-DIR~26kLcxp`+2`kG#=coy6TTUE_reQ_wlp4bK(3SNzCAnec;sbeyjWN zpVGrM|BKK>5_fm&T+mv|>hI`;W`(G{^NI_&1Ldw?@LGY6$hacQ7C!!8EgDFbRA6J? zlE2%G{*mF`3Sazbbt`h-tj;)P=s8#4x6v?ra3|ELpEQwK9}RL)%j?6r^eAZ3^U1UV z&dsR`uttg-TXZI5(Y#amBvZKu-kb{p6J>sJ5;4eS;1;^`@}mCGZdJOBYl6N;P;&a* zb*K&#IvuCLK?Mn^i)k}F5f6oT#nj@>L8vc7BJWx1!r$gXmzKxk2e`dzK;|<~hV5h% zugRGEWa?Xry>kUrNfkyhnX|Ek`m=-Cbirbv_IqKq9>v@MjX=MpVJ(#$pBdwryv<~y z!*wp+|NOoOX_`lm#d+A?Lf#UO6m3BTMsk;exj72vh^Y9@!r*&6l23jYI{}#p&gWpx z2ABZaZ3=oCb**{N>G+95bI=cU^@L796I)`b0@(I69>9*Y*-9Qu| zj)76b@S?i8p|zWr&DGCWT(#cEasQ4RxS*hWg#bz>=e<5&$h5E!`9nzLJ+n0#{m=Pk zC@;XFP;a?CkXmZt@sq{=nVav}VP}7lZ!K(0{t==XRQz&QUgL&ZX5TaKLU0jW^?;M$ z---vAn~n2U2@&ryRwlhyeuf)M2jh~JokwB=<`UYj7VbB3Cw=L9aECRysHyc^bCos4 zE>0jscZDPWYXL`w_1?Ir1adCZ31nr=89Dk_%93=iW;}XCdZWh&a3h9C@r6$yj zbt{>=*sZ!8nA4llzx2&T9Z-ckx5eb zM2S%5z|Q;F_r1Ujp7CCM{{aS#{Q%*ytg7Gf4pst^Frarbs6>H^`Js9k)pr2zK(K8) zWd!1C`nI4xM%EJCmBbvIIqT zAJ|?nwiK`mDzyhSQ@?C4)R|}-N1Qw*1t1#t;?eWwvbxmmS;>+s2Fh%;fr4P&j_|(I zUqQ?9^kqdUl}<)AnX3iPC}5ls1i|Z=uR6}(bWeU8TDwR%0-1xgaffj972sJwGRENk z6qKQdLh#pF(T;hWL-auL0o0xXY@%8`5V5xq$k|Sfg9t@&FDRsP`UH~THT@dz)?$K1 zc|#Wn#B1;|QJvp23es)fzO!n7V~Iq|!%<&iERvoDEzp`6eE`#Kw)JAO*un)2SOD01 z)+5MJ>{5O^7$ptn8-h*Stj!07OkNS18u6PqnzvJG;QHJ>(1*P9(o8H8%L50aLya=f zy!Z1YmPUI=gA2uyjd2M}Ol-J6!=!C{^KWF*@f@j$k4nOr-7#~}UDgL)S6a5C1*#EH zv_SIl$vbj;e;Zbe`fDWrfX;Di;a_GE>$k>Up=d>%#)YKs7ma3C*#saZMjt>CiS`24 z_Ib)LFliHhME8Hlwb{(YFW=8O4H9+LU~=|P#D~Bkf$(|hP9k)9VX6FbcJ#5HaWFd@ zAbKF+o7L~2y%?fVI8Al)*tdP70ZiF^58|%3S^xWp;b9@VEl0gBAI4Ry+%*2YG3OP9 zT&CgJ8@z1f+BN3yR>nR8+VUoDnToy;E}n#}CqIJa*01S5@zDS|<(%)apc}G9{D8DZ!Dz@ysQ$W|QAeFv3 zCKqax2If2i3j+pm0?-egU7*MwsN!i8=eDNeN-&A5R;&e1_$sCbrVG=UH z+GXqoO%eN&M6sW`Acm^qJ%8-z(X;10E4p)^zY`x_YtM~6X888VLT+aH#G(4<&)({t zm$uCdsbg)3%?cgbH3V01r9!EeXA>AFSQN8glABF6zzXXv4X@WkeH@2AxeP~1>58>J zg|OEQ=uQSMplsU67sjgukTA*JPMGMn7dKn&eiLe|wK>TXC>qpZ0a$^_NYM0c5h~i) zzK{eW<4_m_bpYC&#d1eROlO<`7p*ZbZN zyzhJa9yscMs{z|pV-rejkVgMlEGc^9-sh!Nr@KE!sIBo9DKgr&?*rjHKmSiCi|-f4 zqOY;GC@sGT?-6K&K8_yLPJ`#EBY)CPQ;}4kicyO|_GBjrM8-Wi$|`-LG{oiwLxnH} zN2zt!n681b#411c_gWxpD&X9UVVJBb|5+z3@m{b1>lxXxUWX>AU?4qqYVi52&x zitQNc|A^=;D(b>lwD!V&LJE?JlQl4;V0q0bQt)`pqLST)0L+pS#=^y^mCs<6K)QDP zp?RWbs}efG0{RGjB+?b9n<@}4(v{4SO8I!AQD6!WC{(@1%*ge`KmpcDYjk)OwQ2$v zb}@-96^wNLjkuzc!h}ImhG}s&_J*K=XUD`IJz^Dw)o8VU7Tz$sTtw&ZA-O=|*ccI- z98PD%TXOFAblM>|*Tt}gWAa#6fd}G#YD|aV1M5);zIC5->C0ND(OrKY8w$ z-A!kXFVZZYn%}WgLFAO)0Ca?MNtm@my7l=z6|KU7tFn!W>w_L4LkpPQj6cl3O)Y@( zMx#m75xdaB{X>;WJBFsn1ue+*^v;zRL%B*4clQ;lO59~Ka;Vo8icnsW7^GG|cha&w z%vmKi?DyOyi*ga|XplGI0g? zh41HIAvW4xEtR9CXB~EhvA`z-okhMkt9-R{gEw@#iw2IpDs9cnRy*J~!!4z&rJHL~ z*n=t`A&uw$=pE0!x36%>Iu&0_)+gYv$vd9)QOIHGmh334sH1P`pEI$ROII1r(0EfF zH7~^<=Odzy0HP6x{BfmVOeuBlMfkPZk4k+klIjS9bt-nC5fKwUc~D!gJ5wN%{uRpY z=Ow-mcg$*8Mw#OC9$u(g}5 zvoh%fULp%g?fqV&K#U{z|3ni-lg>N#!d{ypR6(UNex*z?0-2qYEf8sllsMU#sMlUv z(N!?LUs?j@E4InO(Ef~01^hqaV8<^U1S*;MkC>(sfyCRR`}K*thc>k;Qn`nv0>_Tn z=w|rhFxj^z+<3MC*v(T;A&Gc>Q)yl5+lw^s{g1B zd7`QfsrtAMkAOGhJ&rH4mt|_lGX74Vr0~x2OS0+g(Ny|l8V~%z79%orE&T8#;f6=dgAbQzh+--K;a#@a&>kDrGhU)jHLqL1Tx}Wyc1Lm|E``E%@?4yxI`TP3 z>Z*xC)7y05W)Ag>UR-1SD)DY8=w4u(ZatOYt%GecAP6=ll78=Dela|9oPnsprBD2g ziP{8hAoT(_L8YN*sTRVvY9=yiq=zyOvmEQK4*8;jbQ z@9uUU)L}U)wrFf>JX>xpUbQ=W*3D??$(5B@AIfPXI?lDPbe{*-04xnkuOoUga7S3fQ411j$@MRV*BFRHkV2W16k1Ky3j_lGApoD#WHoucxmz6E`bk zkZJSK@W8Go7ngiC22~Nk%>%QfnHQGh6D@9EwD(7F+W3Wb@luE6eAKde7S&sC_9X|e ze9E3BYv;m=o(kMT=j5iiJb5&Et?)3^f^G;H&oRdFC4ZsvZf&+5r`)QvZ(a6Q_fK#mxfl zM7=jmMPlMG;%zz{cBQv_f0v@_J|}%wzxoUtc6ajoGAofLi(A-W#IdPRQmHn!nUt&d z?%s{U)E|m0g^htj?4R0z7tiS3;3sImr{5nhKd}4a>ySDxMdonv-p^GEbKHHOy81hXJpFS(jc8l23)0BqM zN#c0|96She(p<5_@p!lsWoNS5X?v);Yqq z{4m!5kOE(CE#v|13(@K7v6|XV10kK~T9gXzDnJJlK=RzVRj2wQ9yQgWDHwi>Gyab$2V4KtmBL^q@Q%ExE3B|*&8XPvMN6GSc z32x#7Tbv+3kMW;14P2MlYe&iKTT)V6Pm7K~V)CW3W(?tA0k(}zGGCJWdEo>>;P1xC zlGuRr-Xa+pVo`@Z3R5RvPr}5@qJ?6<;(U%2alZq@kW- zvX-5v5R8j`A`C24zEHVx)TVOR>n{a&%_E1-O#lm>ENI&PFq5Yv?rB6hAOFc8gZ?U>5FzvbHD zFac2oN$*VCV484^KEXizyg=mI(o793;e8`n*fXTQ(>bs3-K1VukZ=MDuFqN2zb9;_ zE)aE5E?wY}2v4(SYD=2PC`MuOIXFLUQ7dpi-6J7uOeu8uu^$AZ0k1mTm!RqXg`nFr8$sR5gA3nLgr#~o@h;gsU;#yDqc0!_#$yx7Tr0k~5`(MZt~{*O9o_U8PShaZ zbc$s+Y>rf(AF15?CoDmPz?LOQ|2s<%LH7?V0q_imy;mmQ8dYpiSB)x`zr`9=tW2gq zrQcW@!N!MC>BB~JRPjsGr89}#iLEZo1TfX{eV0pxj;G5Zf04 zT&qgCoO@1asciuIuh|&t@>pBri|ASV$n#U56?hSzOBx+-yEGrA9sboXmdwtI2R4($ zZAlL&6NLWpDq5m_Sv~)PeYpGJi@)lsz`+V(m(b{xYv9l-a1}Vqe{2TLuU}+D`~}A4 zZ1qQ>NYIvHN&Se?EeFd-!kz`pAV^e2(8vD6ZRY z9&dCX@K==@#Q+!Y6rtSlKbZAmk`_VXFL0JBgWeLm;|1m60uS_!{VLl!t7gRb0^<2^ zDn8=1l`*LrwS5wOjD(CpdV4yW#J2Z-bw7+~FtsYx)zc+sPF{{`&*jMq)x+$MB>mJ0vRD*I$Cp4FHifI0BRafbH|Dbu;zr zvP0HZdKWwPcx%b-hOOB(BX^!sm6&+fTkoaA^JoaT(o@7noKCNTtxzX81v_mKobwJ$ z2!G+7HS*EU4ImlqySoZq@XSR#={<|xKliv~x!b+$3FNV%&HR*Sq$~I3h_%07HUm*x z92|IJQY>!ao_Q7LtnaazSxK2PcOUf_>E~76Tt=7qoKouh!E5LU?30qXJSFm2#*=}B zL@*o8hdM)jtIS%`^SC6LKdIq!*UZM!%kKb_A)E@1Be-~UdsPJsvQy`TLbZJ1YZ;CQ zY?g)sBkA_qJ5US3!dSR1f6>9oaI+5CX@yr|@wH92f1EgWO4=GRnL!mR z6Dth)+~)Lx<5V^PJ;Iq)l3D&1vBuIm{F7RY#$aD%04}t!?X6GvaL){wSnV~(lYjN z>9Zr{rm0|fv)=z zI^KaTmHHxT-s4H9D-;@sCx=BSOH>I({H!GZfRVWYkDoo?Br-L_70Nar{iLRip3P@1 z>>_P7OVdr&mt#l&%TCW7}*2&KGUifz3NOsVqGHT>?W&Q5)9-oe!`a%YqY&BfW z&^$P+WOV=m%Hp(zYZXjn@)gYJlLuuBnBfHa+pjI+%{oTe^H`1)8r`2V;4SD=&i2Ig zuKL+`XDALoB+!5qK+Y}B7XGAf4@3%(aP>*|)}OlK^VGeu5l695gc{cg!v6GV6D`x( z{hy<0j(Q~)co~F6UTMjCeDH>3RWz@yG_nfFlRs%=C{AHnH3WMzatM62c^FM*+?@ zNCJ{?Al2)hVf9(>vem%*$!aGnPWK$qmJ7D;P7$ItOdUXws#Q=yE|AxBZUUb_fi z06a@XUF>>A+ z&5O&xsUSTcFpU(bht)*8TPsbS*M=A0KwoAlystI>MJUEmrr~>Sz?8PJ2VM~^^;u+q z^9R-V#W_;=`>TN5@Ampg^jYka=X0Y~F|J0~YATQzOYng?11367> zzIWuuYa4+;=-edb@vHyr4ZX-M7X(t3hQ@iHsD#uWQQyLtl>80EQKS3gx%GuX^V+HZ zl4VG+Cv6$Hf3M*F80^s-Q&>f%=rj08(ca0yykf#cEMu&S3ni;4NB{3Glii}lU4Rm{ znNv7l(0)PP>Qu3MY$tzG)|jaJx{S!^j}xhL&ww>{SLWWz7yVfaX2;T_XK^8a#Cq4# z`NLCVIMA-O%YA($J%ynhr1KXB@_cVrwl~e?4fV9r5dT&6Bb}f<*0~Fv&u@ zxSLs|UI+|A@T%FMH|SkBJ3w8H1K;gA-@m+x6^QJ%m)hyH##IQ?#=aDa#zqi#nM}h# z;h?N*`p3Z`-chH6gM?O^)qj=j!GaR2<1kTW*U#ZZRKGSajonChvm^ITP1{NQe0m^I^3&(JlacIYZLOA|zWbCOWC=nu~> zuyJ5<{HS|O6JxlaM&?}2lzDn_FSjC2~xgT6F@Nh(; zl;q^L9-G%OyTIcYQVZ+wAtlaqBL~RHHl`1CZWkm3u~!>L<~}^*y1cTy{TTFqb~CQR zaH^h=SJ91?ELC^&i-N5q;o@$h5E^I#rylmZT)}fA8jI(_x$=rMsRTe9&D0eMCfxtxlN`{Xv(@-*?}M^4NL?P{5duLX10E2*qCXnXuUDHsiY zEP9D%Lzfe7-f+raFLIfkGU(Wthyf9n(aJH>qh%0LG2>35uVCRttj%|68!!7(sna#Ta>%1OgYMptv|btYR05F07~|r*=id-m_-9uhb=; zv^W`R*WI949m?W?qtChHrE0@C)Lb9RVrxUASQjK&ALt<|={zJRIk$4%{F}^A3v*(& z19whY-FT{8bMdDolPY4H^W{Myzeo4)b*95Ej6b`Q6101fy1x1p^pCHT!n1%_ihB3t zuuzkJVNVcitD;`ev5eEM1lvmdYGvkC)^5p8OyvS?1?ix*LHg(J8=Bccr7*-ETTe(g z#Yn4cMAjO&O}+qmul-_L3yaHd@pO=nc{Q41@qRI za(CS>#7J>sF3uJvsYJtGn+vvSt+LiETi$=8s|uJesTgEtXqu$;7QtBef^kb8mcDJ) z85yD2z{f0*lWJaxBe2#yrJZ)N5OYD%+Ib20F3mFK4xIqouRq*Cn&nJ@7Al?qE%b$V zz2dRjs&1*ZS1J#k{f>=1rqRJ1qTEN_%)9DUf5KrRrwie4x6dvf%ky34#E=WwJOy4{ zR3eEGK~Gh0{p^*~h+Zh?Px3Im%xJW)82_?>Ys1-e>xbV=R^fzp(k0_~@faW7JY}}Un#Nn(`i35%qdRbrQ+`Q;BXd+FcSo-<;#GxJ8;5Cadq$YSc*3yLQfRxwm2V@@zF6dZ ztUo9e$(nM;?%A5;F`b%yQnwOR9?07FFu#2;e8S-ow){1hoX1CV#AU{|%X_uRsN5;^ zX}L}m_a%;O-94iRFJ)z-@UT5reGjbx$S8z$j#n>H!>ifH&?W ztk*9Sb>x)VZ*2Z9sq=04Q6zF|3U`8L80Gow&=UJ8x>!oC5nZ>!Lck*$pA@1*ehSu? z=*#U{I8uCa?)<9t!ToX~VZi3ORn@|#9Wq3w2CUFKq{C3G-ARSf7t0ii5kyXl3^aP` zuAE78zQd3O${g(}kX6r+UryeUuV_P)iP0YRLvhNxp6V^=bsRrvM?gZX89qC=EHCK!W6Lx8` zJrJ;!P3|^1cjonM>2sZ~z{B#ch*{@^ev%dX?qGFUmI0$x2q~}I#H0n)S#l0AYR?VN zVsIoV;7M-}W+_03_|LxhSjPY)vV5dAuP)GIILx59E)-ww1Jtt2fAvjT`1=&MmA4vn z=jggAc}4GOQfkSU8BJJv@N_*#9hft2p~Z%?+0Nu4iaxc@DO>X!v3Lu@{8^`~Ki;0M zqqo!Ww4~<23~}70cDvie*SP9QgW=1?3g3p0yNz=EX%F|?7axH8Kfk%cYH&T{S6Aq7 z9BHQhKRMKz>b13nwn`9rXNmXT3*KUg^aS5O{S@Bqvjk-k5bJC;XvhyEbCfowT`KeH zj7>IHzzlOPAJFXYG-wU9rZFOQ^C>A8@`9=bfUFA#fKokT=JTJ`_XLF(rt2G%>Dh+e z9%aVQlIkC^h8>d>^GABG&EcK~twH{Ep{;NhKi>eZ(+$uuj1xipmg_9561&ZwoBQOt zEHXLe8HdD+ogfM3^B=p9_sLilqtE|(fC|4Je2|aplKE}M^S)POr5KoaW{I^vUbVBI zRprX*Sh;YOoQymtn(^G=1-;)ncn?5tl3~x>ph3RwzrM4qc)b2Tn~=QpgOB;F&NRW! zn7mK@rxw}PDfBxq0!K4NIRjbcrKLO?f`F-%cioqjxv_`ZDeQk_-^0Dizx!o9G}~YnJW1H?HYz@Cigo=hpccosL)DqZkpxR9Xe3w!?(18`Al-mr63*P`(^*ii`9T|q-ES;fOR?D4KNYDwS%2xeU_tSx;VX=aKW`B{ zU2agNumiVW`p}e*!(b8fRHi(PiO=88mJ+hH<6QcB8!?nd0Z-wAjtJzGkAJgJWKHBo%0RBehy(SDfmOo&7X?7&+x-YzsCX}MZx3$SS%|9M*E#N{(w@)3jz|*5C!S8Of zH1a3^BC)iWj#UZLRD*k{9zYI%UIJTXBq&Leqc*6gkFAv_vD$(F^#4^2e^cgvf@^>K l;=iKxzr97-dwFz*#wnWZ*&24hlm`4kUcGT8OWydwe*vPLE(ibs literal 0 HcmV?d00001 From 20b9e46ab7a1e021437dcc6fe807e2095db199b3 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 28 Nov 2023 15:02:25 +0000 Subject: [PATCH 25/37] Add release notes --- .../v6.9.0/Diffraction/Single_Crystal/New_features/36481.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36481.rst diff --git a/docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36481.rst b/docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36481.rst new file mode 100644 index 000000000000..8f7aa9bfb38f --- /dev/null +++ b/docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36481.rst @@ -0,0 +1 @@ +- New shoebox integration algorithm for single-crystal Bragg peaks :ref:`IntegratePeaksShoeboxTOF ` \ No newline at end of file From c04b1e709ecb76543e547a18c8dbd1120e2b713e Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 28 Nov 2023 16:39:15 +0000 Subject: [PATCH 26/37] Always use QLab angle for closest peak in all cases & related bug fix Also scale TOF extent of shoebox by ratio of peak TOF if GetNBinsFromBackToBackParams=False --- .../algorithms/IntegratePeaksShoeboxTOF.py | 32 ++++++++----------- .../IntegratePeaksShoeboxTOF-v1.rst | 11 +++---- 2 files changed, 18 insertions(+), 25 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index ac2e09707ef1..ad9f291bcbbd 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -384,6 +384,7 @@ def PyExec(self): 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) + fwhm = fwhm if get_nbins_from_b2bexp_params else None # not calculated but not going to be used weak_peaks_list.append(WeakPeak(ipk, ispecs[ipos[0], ipos[1]], x[ipos[-1]], fwhm, ipks_near)) else: if status == PEAK_STATUS.STRONG: @@ -400,11 +401,6 @@ def PyExec(self): set_peak_intensity(peak, intens, sigma, do_lorz_cor) if len(ipks_strong): - # set function for calculating distance metric between peaks - # if know back to back params then can scale TOF extent by ratio of FWHM - # otherwise just look for peak closest in QLab - calc_dist_func = calc_angle_between_peaks if get_nbins_from_b2bexp_params else calc_dQsq_between_peaks - for weak_pk in weak_peaks_list: # get peak ipk = weak_pk.ipk @@ -415,26 +411,30 @@ def PyExec(self): ipks_near_strong = [] for ipk_near in weak_pk.ipks_near: if results[ipk_near].status == PEAK_STATUS.STRONG: - ipks_near_strong.append(ipk_near) + ipks_near_strong.append(int(ipk_near)) if not ipks_near_strong: # no peaks in detector window at any TOF, look in all table ipks_near_strong = ipks_strong # loop over strong peaks and find nearest peak ipk_strong = None - dist_min = np.inf + angle_min = np.inf for ipk_near in ipks_near_strong: - dist = calc_dist_func(peak, peaks.getPeak(ipk_near)) - if dist < dist_min: + angle = calc_angle_between_peaks(peak, peaks.getPeak(ipk_near)) + if angle < angle_min: + angle_min = angle ipk_strong = ipk_near strong_pk = peaks.getPeak(ipk_strong) - # get peak shape and make kernel nrows, ncols, nbins = results[ipk_strong].peak_shape + # scale TOF extent if get_nbins_from_b2bexp_params: # scale TOF extent by ratio of fwhm strong_pk_fwhm = get_fwhm_from_back_to_back_params(strong_pk, ws, strong_pk.getDetectorID()) - nbins = max(3, int(nbins * (weak_pk.tof_fwhm / strong_pk_fwhm))) - nbins = round_up_to_odd_number(nbins) + ratio = weak_pk.tof_fwhm / strong_pk_fwhm + else: + # scale assuming dTOF/TOF = const + ratio = pk_tof / strong_pk.getTOF() + nbins = max(3, round_up_to_odd_number(int(nbins * ratio))) kernel = make_kernel(nrows, ncols, nbins) # get data array in peak region (keep same window size, nshoebox, for plotting) peak_data = array_converter.get_peak_data( @@ -467,7 +467,7 @@ def PyExec(self): # plot output if output_file: prog_reporter.resetNumSteps(int(len(results) - np.sum(results is None)), start=0.0, end=1.0) - plot_integration_reuslts(output_file, results, prog_reporter) + plot_integration_results(output_file, results, prog_reporter) # assign output self.setProperty("OutputWorkspace", peaks) @@ -490,7 +490,7 @@ def round_up_to_odd_number(number): return number -def plot_integration_reuslts(output_file, results, prog_reporter): +def plot_integration_results(output_file, results, prog_reporter): # import inside this function as not allowed to import at point algorithms are registered from matplotlib.pyplot import subplots, close from matplotlib.patches import Rectangle @@ -518,10 +518,6 @@ def calc_angle_between_peaks(pk1, pk2): return abs(pk1.getQLabFrame().angle(pk2.getQLabFrame())) -def calc_dQsq_between_peaks(pk1, pk2): - return (pk1.getQLabFrame() - pk2.getQLabFrame()).norm2() - - def get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox): x, y, esq, ispecs = peak_data.get_data_arrays() # 3d arrays [rows x cols x tof] x = x[peak_data.irow, peak_data.icol, :] # take x at peak centre, should be same for all detectors diff --git a/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst index ca32e8a34ee0..fc7048d905c5 100644 --- a/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst +++ b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst @@ -41,13 +41,10 @@ The algorithm proceeds as follows: shoebox dimensions from the nearest strong peak (``WeakPeakStrategy="NearestStrongPeak"``) When looking for the nearest strong peaks for ``WeakPeakStrategy="NearestStrongPeak"``, the algorithm first checks for -peaks in detector IDs in the data window around the peak, if one isn't found the nearest peak search depends on the -parameter ``GetNBinsFromBackToBackParams``. - -1. If ``GetNBinsFromBackToBackParams=True`` then the closest peak is defined as the one with the smallest angle between - QLab vectors and the TOF extent of the shoebox is scaled by the ratio of the FWHM of the weak and strong peak. - -2. If ``GetNBinsFromBackToBackParams=False`` then then the closest peak is defined as the one nearest in QLab. +peaks in detector IDs in the data window around the peak, before looking at the whole peak table. +The closest peak is defined as the one with the smallest angle between the QLab vectors of the two peaks. The TOF extent +of the shoebox is scaled by the ratio of the FWHM of the weak and strong peak if ``GetNBinsFromBackToBackParams=True``, +otherwise it is scaled by the ratio of TOF (i.e. assumes dTOF/TOF resolution is the same for both peaks). Optionally if ``OutputFile`` is provided a pdf can be output that shows the shoebox kernel and the data integrated along each dimension like so From cb21c005b0b525e26ce81e1fc36b2bb6a3f648e5 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 6 Dec 2023 09:19:23 +0000 Subject: [PATCH 27/37] Update doc test value --- docs/source/algorithms/FindSXPeaksConvolve-v1.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/algorithms/FindSXPeaksConvolve-v1.rst b/docs/source/algorithms/FindSXPeaksConvolve-v1.rst index 34c86471d329..e061492541e3 100644 --- a/docs/source/algorithms/FindSXPeaksConvolve-v1.rst +++ b/docs/source/algorithms/FindSXPeaksConvolve-v1.rst @@ -54,7 +54,7 @@ Usage .. testoutput:: exampleFindSXPeaksConvolve - Found 323 peaks + Found 307 peaks .. categories:: From 452cde24ea2e8571a6ad73324d304dabc1b65cc9 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 15 Dec 2023 16:04:11 +0000 Subject: [PATCH 28/37] Add refernces and background to documentation --- .../IntegratePeaksShoeboxTOF-v1.rst | 39 +++++++++++++++---- 1 file changed, 32 insertions(+), 7 deletions(-) diff --git a/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst index fc7048d905c5..0dc873d3fc33 100644 --- a/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst +++ b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst @@ -6,16 +6,35 @@ .. properties:: +Background +----------- + +Fixed mask/shoebox integration methods have previously been developed for Time-of-flight (TOF) neutron diffraction [1]_. +For monochromatic neutron diffractometers, a dynamic shoebox method that optimised the shoebox dimensions on a 2D +detector by maximising the Intensity/sigma ratio of each peak has been shown to improve the R-factors obtained [2]_. + +This algorithm generalises the 2D dynamic shoebox integration of [2]_ to 3D peaks in TOF Laue, +and several improvements have been made to suit TOF Laue data + +1. Avoids nearby peaks that are closer to another peak position when optimising shoebox position + +2. The shoebox dimensions of weak peaks can be obtained from nearby stronger peaks with scaling of the TOF dimension +to account for the instrument resolution) + + Description ----------- -This is an algorithm to integrate single-crystal Bragg peaks in a :ref:`MatrixWorkspace ` by -summing up counts in a shoebox of size ``NRows`` and ``NCols`` on the detector and ``NBins`` along the Time-of-flight -(TOF) direction which can be specified in one of two ways: +This is an algorithm to integrate single-crystal Bragg peaks in a :ref:`MatrixWorkspace `. +An initial integration is performed by summing up the counts in a fixed shoebox of user specified size ``NRows`` and +``NCols`` on the detector and ``NBins`` along the TOF direction. The dimensions of the shoebox are then optionally +optimised to maximise I/sigma. + +The default ``Nbins`` can be specified in one of two ways: -1. Provide ``NBins`` - number of TOF bins in the kernel +1. Provide ``NBins`` directly - number of TOF bins in the kernel -2. Setting ``GetNBinsFromBackToBackParams=True`` and providing ``NFWHM`` - in which case ``NBins`` will be NFWHM x FWHM +2. Setting ``GetNBinsFromBackToBackParams=True`` and providing ``NFWHM`` - in which case ``NBins`` will be NFWHM x FWHM of a :ref:`BackToBackExponential ` peak at the center of each detector panel/bank at the middle of the spectrum. @@ -27,14 +46,14 @@ the background shell as in the peak region. The algorithm proceeds as follows: -1. Take a window of the data, ``NShoeboxInWindow`` the size of the shoebox (``NRows`` x ``NCols`` x ``NBins``) +1. Take a window of the data, ``NShoeboxInWindow`` the size of the initial user specified shoebox (``NRows`` x ``NCols`` x ``NBins``) 2. Convolve the shoebox kernel with the data in the window to find the position with largest Intensity/sigma (closest to the predicted peak position than to any other peaks in the table). 3. Integrate using the shoebox at the optimal position -4. If the peak is strong (Intensity/sigma < ``WeakPeakThreshold``) and ``OptimiseShoebox=True`` then optimise the shoebox +4. If the peak is strong (Intensity/sigma > ``WeakPeakThreshold``) and ``OptimiseShoebox=True`` then optimise the shoebox dimensions to maximise Intensity/sigma 5. If the peak is weak, it can be integrated using the initial shoebox (``WeakPeakStrategy="Fix"``) or using the @@ -80,6 +99,12 @@ Usage I/sigma = 100.49 +References +---------- + +.. [1] Wilkinson, C., and A. J. Schultz. (1989) J. Appl. Cryst. 22.2, 110-114. + +.. [2] Wilkinson, C., Khamis, H. W., Stansfield, R. F. D. & McIntyre, G. J. (1988). J. Appl. Cryst. 21, 471-478. .. categories:: From ee1c5895ad20528967804163e6c55c6713c2f041 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 18 Dec 2023 18:25:22 +0000 Subject: [PATCH 29/37] Account for bin width when scaling TOF nbins for weak peaks Also defined WeakPeak as a data class --- .../algorithms/IntegratePeaksShoeboxTOF.py | 42 ++++++++++++------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index ad9f291bcbbd..42572ace34f6 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -22,6 +22,7 @@ EnabledWhenProperty, PropertyCriterion, ) +from dataclasses import dataclass import numpy as np from scipy.ndimage import uniform_filter from scipy.signal import convolve @@ -38,13 +39,14 @@ class PEAK_STATUS(Enum): NO_PEAK = "No peak detected." +@dataclass class WeakPeak: - def __init__(self, ipk, ispec, tof, tof_fwhm, ipks_near): - self.ipk = ipk - self.ispec = ispec # spectrum at center of kernel that corresponds to max I/sigma from convolution - self.tof = tof # TOF at center of kernel that corresponds to max I/sigma from convolution - self.tof_fwhm = tof_fwhm - self.ipks_near = ipks_near + ipk: int + ispec: int # spectrum at center of kernel that corresponds to max I/sigma from convolution + tof: float # TOF at center of kernel that corresponds to max I/sigma from convolution + tof_fwhm: float + tof_bin_width: float + ipks_near: np.ndarray class ShoeboxResult: @@ -347,8 +349,7 @@ def PyExec(self): nrows = self.getProperty("NRows").value # get these inside loop as overwritten if shoebox optimised ncols = self.getProperty("NCols").value ispec = ws.getIndicesFromDetectorIDs([detid])[0] - itof = ws.yIndexOfX(pk_tof, ispec) - bin_width = np.diff(ws.readX(ispec)[itof : itof + 2])[0] # used later to scale intensity + bin_width = get_bin_width_at_tof(ws, ispec, pk_tof) # used later to scale intensity if get_nbins_from_b2bexp_params: fwhm = get_fwhm_from_back_to_back_params(peak, ws, detid) nbins = max(3, int(nfwhm * fwhm / bin_width)) if fwhm is not None else self.getProperty("NBins").value @@ -361,7 +362,7 @@ def PyExec(self): peak_data = array_converter.get_peak_data( peak, detid, bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge ) - x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox) + x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, kernel, nshoebox) # perform initial integration intens_over_sig = convolve_shoebox(y, esq, kernel) @@ -385,7 +386,7 @@ def PyExec(self): # 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) fwhm = fwhm if get_nbins_from_b2bexp_params else None # not calculated but not going to be used - weak_peaks_list.append(WeakPeak(ipk, ispecs[ipos[0], ipos[1]], x[ipos[-1]], fwhm, ipks_near)) + weak_peaks_list.append(WeakPeak(ipk, ispecs[ipos[0], ipos[1]], x[ipos[-1]], fwhm, bin_width, ipks_near)) else: if status == PEAK_STATUS.STRONG: ipks_strong.append(ipk) @@ -432,22 +433,25 @@ def PyExec(self): strong_pk_fwhm = get_fwhm_from_back_to_back_params(strong_pk, ws, strong_pk.getDetectorID()) ratio = weak_pk.tof_fwhm / strong_pk_fwhm else: - # scale assuming dTOF/TOF = const + # scale assuming resolution dTOF/TOF = const ratio = pk_tof / strong_pk.getTOF() + # scale ratio by bin widths at the two TOFs (can be different if log-binning) + ispec_strong = ws.getIndicesFromDetectorIDs([strong_pk.getDetectorID()])[0] + ratio = ratio * get_bin_width_at_tof(ws, ispec_strong, strong_pk.getTOF()) / weak_pk.tof_bin_width nbins = max(3, round_up_to_odd_number(int(nbins * ratio))) kernel = make_kernel(nrows, ncols, nbins) # get data array in peak region (keep same window size, nshoebox, for plotting) peak_data = array_converter.get_peak_data( peak, peak.getDetectorID(), bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge ) - x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox) + x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, kernel, nshoebox) # 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) # scale summed intensity by bin width to get integrated area - intens = intens * bin_width - sigma = sigma * bin_width + 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 @@ -490,6 +494,11 @@ def round_up_to_odd_number(number): return number +def get_bin_width_at_tof(ws, ispec, tof): + itof = ws.yIndexOfX(tof, ispec) + return ws.readX(ispec)[itof + 1] - ws.readX(ispec)[itof] + + def plot_integration_results(output_file, results, prog_reporter): # import inside this function as not allowed to import at point algorithms are registered from matplotlib.pyplot import subplots, close @@ -518,11 +527,12 @@ def calc_angle_between_peaks(pk1, pk2): return abs(pk1.getQLabFrame().angle(pk2.getQLabFrame())) -def get_and_clip_data_arrays(ws, peak_data, pk_tof, ispec, kernel, nshoebox): +def get_and_clip_data_arrays(ws, peak_data, pk_tof, kernel, nshoebox): x, y, esq, ispecs = peak_data.get_data_arrays() # 3d arrays [rows x cols x tof] x = x[peak_data.irow, peak_data.icol, :] # take x at peak centre, should be same for all detectors + ispec = ispecs[peak_data.irow, peak_data.icol] # crop data array to TOF region of peak using shoebox dimension - itof = ws.yIndexOfX(pk_tof, ispec) # need index in y now (note x values are points even if were edges) + itof = ws.yIndexOfX(pk_tof, int(ispec)) # need index in y now (note x values are points even if were edges) tof_slice = slice( int(np.clip(itof - nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), int(np.clip(itof + nshoebox * kernel.shape[-1] // 2, a_min=0, a_max=len(x))), From 77517afa8492b1db1e783f52f9e20c41f354ef19 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 18 Dec 2023 18:26:42 +0000 Subject: [PATCH 30/37] Break out of sorted loop when below threshold --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 42572ace34f6..7a78e5d9a1c5 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -585,6 +585,8 @@ def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol else: continue # executed if inner loop did not break (i.e. pixel closer to nearby peak than this peak) break # execute if inner loop did break and else branch ignored (i.e. found bin closest to this peak) + else: + break return imax_nearest # could be None if no peak found else: # no nearby peaks - return position of maximum in data From a3e01e9bb812f1c2910ecf74608e511a30ce2bc2 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 18 Dec 2023 18:28:02 +0000 Subject: [PATCH 31/37] Rename cen variable in plotting code --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 7a78e5d9a1c5..a123c2dc63b7 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -111,16 +111,16 @@ def plot_integrated_peak(self, fig, axes, norm_func, rect_func): fig.suptitle(self.title) def plot_shoebox(self, iax, ax, rect_func): - cen = self.ipos.copy() - cen.pop(iax) + center_peak_pos = self.ipos.copy() + center_peak_pos.pop(iax) # plot peak region shape = self.peak_shape.copy() shape.pop(iax) - self._plot_rect(ax, cen, shape, rect_func) + self._plot_rect(ax, center_peak_pos, shape, rect_func) # plot kernel (incl. bg shell) shape = self.kernel_shape.copy() shape.pop(iax) - self._plot_rect(ax, cen, shape, rect_func, ls="--") + self._plot_rect(ax, center_peak_pos, shape, rect_func, ls="--") def _plot_rect(self, ax, cen, shape, rect_func, ls="-"): bottom_left = [cen[1] - shape[1] // 2 - 0.5, cen[0] - shape[0] // 2 - 0.5] From 674f4388125e46fe3324c7bcf0a5867285fd5df4 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 19 Dec 2023 19:33:14 +0000 Subject: [PATCH 32/37] Refactor PyExec to reduce length --- .../algorithms/IntegratePeaksShoeboxTOF.py | 93 ++++++++++++------- 1 file changed, 60 insertions(+), 33 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index a123c2dc63b7..5fc43cb73168 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -335,7 +335,6 @@ def PyExec(self): prog_reporter = Progress(self, start=0.0, end=1.0, nreports=peaks.getNumberPeaks()) for ipk, peak in enumerate(peaks): prog_reporter.report("Integrating") - status = PEAK_STATUS.NO_PEAK intens, sigma = 0.0, 0.0 detid = peak.getDetectorID() @@ -363,24 +362,27 @@ def PyExec(self): peak, detid, bank_name, nshoebox * kernel.shape[0], nshoebox * kernel.shape[1], nrows_edge, ncols_edge ) x, y, esq, ispecs = get_and_clip_data_arrays(ws, peak_data, pk_tof, kernel, nshoebox) - - # perform initial integration - intens_over_sig = convolve_shoebox(y, esq, kernel) - - # identify best shoebox position near peak ix = np.argmin(abs(x - pk_tof)) - ipos = find_nearest_peak_in_data_window(intens_over_sig, ispecs, x, ws, peaks, ipk, peak_data.irow, peak_data.icol, ix) - - # perform final integration if required + ipos_predicted = [peak_data.irow, peak_data.icol, ix] det_edges = peak_data.det_edges if not integrate_on_edge else None - 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) - 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) + + intens, sigma, status, ispecs, ipos, nrows, ncols, nbins = integrate_peak( + ws, + peaks, + ipk, + kernel, + nrows, + ncols, + nbins, + x, + y, + esq, + ispecs, + ipos_predicted, + det_edges, + 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) @@ -409,22 +411,7 @@ def PyExec(self): bank_name = peaks.column("BankName")[ipk] pk_tof = peak.getTOF() # find nearest strong peak to get shoebox dimensions from - ipks_near_strong = [] - for ipk_near in weak_pk.ipks_near: - if results[ipk_near].status == PEAK_STATUS.STRONG: - ipks_near_strong.append(int(ipk_near)) - if not ipks_near_strong: - # no peaks in detector window at any TOF, look in all table - ipks_near_strong = ipks_strong - # loop over strong peaks and find nearest peak - ipk_strong = None - angle_min = np.inf - for ipk_near in ipks_near_strong: - angle = calc_angle_between_peaks(peak, peaks.getPeak(ipk_near)) - if angle < angle_min: - angle_min = angle - ipk_strong = ipk_near - strong_pk = peaks.getPeak(ipk_strong) + ipk_strong, strong_pk = get_nearest_strong_peak(peaks, peak, results, weak_pk.ipks_near, ipks_strong) # get peak shape and make kernel nrows, ncols, nbins = results[ipk_strong].peak_shape # scale TOF extent @@ -488,6 +475,46 @@ def exec_child_alg(self, alg_name, **kwargs): return None +def get_nearest_strong_peak(peaks, peak, results, ipks_near, ipks_strong): + ipks_near_strong = [] + for ipk_near in ipks_near: + if results[ipk_near].status == PEAK_STATUS.STRONG: + ipks_near_strong.append(int(ipk_near)) + if not ipks_near_strong: + # no peaks in detector window at any TOF, look in all table + ipks_near_strong = ipks_strong + # loop over strong peaks and find nearest peak + ipk_strong = None + angle_min = np.inf + for ipk_near in ipks_near_strong: + angle = calc_angle_between_peaks(peak, peaks.getPeak(ipk_near)) + if angle < angle_min: + angle_min = angle + ipk_strong = ipk_near + return ipk_strong, peaks.getPeak(ipk_strong) + + +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) + + # 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 + 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) + 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, ispecs, ipos, nrows, ncols, nbins + + def round_up_to_odd_number(number): if not number % 2: number += 1 From 1dc1e1010e41922087aa8102a44eca4abfdae796 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 19 Dec 2023 19:34:55 +0000 Subject: [PATCH 33/37] Fix bug relating to specturm index of nearest peaks --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 5fc43cb73168..2fcb2c8f34a4 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -478,7 +478,7 @@ def exec_child_alg(self, alg_name, **kwargs): def get_nearest_strong_peak(peaks, peak, results, ipks_near, ipks_strong): ipks_near_strong = [] for ipk_near in ipks_near: - if results[ipk_near].status == PEAK_STATUS.STRONG: + if results[ipk_near] and results[ipk_near].status == PEAK_STATUS.STRONG: ipks_near_strong.append(int(ipk_near)) if not ipks_near_strong: # no peaks in detector window at any TOF, look in all table @@ -621,14 +621,14 @@ def find_nearest_peak_in_data_window(data, ispecs, x, ws, peaks, ipk, irow, icol def find_ipks_in_window(ws, peaks, ispecs, ipk, tof_min=None, tof_max=None): - ispecs_peaks = ws.getIndicesFromDetectorIDs([int(p.getDetectorID()) for p in peaks]) + ispecs_peaks = np.asarray(ws.getIndicesFromDetectorIDs([int(p.getDetectorID()) for p in peaks])) ipks_near = np.isin(ispecs_peaks, ispecs) if tof_min and tof_max: tofs = peaks.column("TOF") ipks_near = np.logical_and.reduce((tofs >= tof_min, tofs <= tof_max, ipks_near)) ipks_near = np.flatnonzero(ipks_near) # convert to index ipks_near = np.delete(ipks_near, np.where(ipks_near == ipk)) # remove the peak of interest - return ipks_near, ispecs_peaks + return ipks_near, ispecs_peaks[ipks_near] def integrate_shoebox_at_pos(y, esq, kernel, ipos, weak_peak_threshold, det_edges=None): From 5a1fab1be9039326ee8fb0f3b72ff32b0a6c9dc3 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 20 Dec 2023 09:55:22 +0000 Subject: [PATCH 34/37] Add shoebox integration to single-crystal reduction classes --- scripts/Diffraction/single_crystal/base_sx.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/scripts/Diffraction/single_crystal/base_sx.py b/scripts/Diffraction/single_crystal/base_sx.py index eef348c0462d..df83d54fffc7 100644 --- a/scripts/Diffraction/single_crystal/base_sx.py +++ b/scripts/Diffraction/single_crystal/base_sx.py @@ -27,6 +27,7 @@ class INTEGRATION_TYPE(Enum): MD = "int_MD" MD_OPTIMAL_RADIUS = "int_MD_opt" SKEW = "int_skew" + SHOEBOX = "int_shoebox" class BaseSX(ABC): @@ -341,6 +342,11 @@ def integrate_peaks_skew(ws, peaks, out_peaks, **kwargs): peaks_int = mantid.IntegratePeaksSkew(InputWorkspace=ws, PeaksWorkspace=peaks, OutputWorkspace=out_peaks, **kwargs) return peaks_int + @staticmethod + def integrate_peaks_shoebox(ws, peaks, out_peaks, **kwargs): + peaks_int = mantid.IntegratePeaksShoeboxTOF(InputWorkspace=ws, PeaksWorkspace=peaks, OutputWorkspace=out_peaks, **kwargs) + return peaks_int + @default_apply_to_all_runs def integrate_data(self, integration_type, peak_type, tol=0, run=None, **kwargs): pk_table = self.get_peaks(run, peak_type) @@ -357,6 +363,8 @@ def integrate_data(self, integration_type, peak_type, tol=0, run=None, **kwargs) if ws is not None: kwargs["ws"] = ws # use this workspace to avoid loading in empty self.integrate_peaks_MD_optimal_radius(ws_md, pk_table, peak_int_name, **kwargs) + elif integration_type == INTEGRATION_TYPE.SHOEBOX: + BaseSX.integrate_peaks_shoebox(ws, pk_table, peak_int_name, **kwargs) else: BaseSX.integrate_peaks_skew(ws, pk_table, peak_int_name, **kwargs) # store result From 2aa6c98d02a301efad2f96e8f13a3525f480cc04 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 20 Dec 2023 11:50:49 +0000 Subject: [PATCH 35/37] Add system test for single-crytal reduction integration --- .../tests/framework/SXDAnalysis.py | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/Testing/SystemTests/tests/framework/SXDAnalysis.py b/Testing/SystemTests/tests/framework/SXDAnalysis.py index a50b3bcf0b7c..7c9afd239c91 100644 --- a/Testing/SystemTests/tests/framework/SXDAnalysis.py +++ b/Testing/SystemTests/tests/framework/SXDAnalysis.py @@ -134,3 +134,39 @@ def validate(self): self.tolerance = 1e-8 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)) From 8d37c2f32f54104aa40ed9d9bb9dfc4ddee0fffd Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 20 Dec 2023 11:53:50 +0000 Subject: [PATCH 36/37] Remove unmatched bracket --- docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst index 0dc873d3fc33..2a7ce4b8a1ab 100644 --- a/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst +++ b/docs/source/algorithms/IntegratePeaksShoeboxTOF-v1.rst @@ -19,7 +19,7 @@ and several improvements have been made to suit TOF Laue data 1. Avoids nearby peaks that are closer to another peak position when optimising shoebox position 2. The shoebox dimensions of weak peaks can be obtained from nearby stronger peaks with scaling of the TOF dimension -to account for the instrument resolution) +to account for the instrument resolution. Description From 76076862515eeaaf28d9b4011b86ac973c9c3dc1 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 20 Dec 2023 12:07:31 +0000 Subject: [PATCH 37/37] Add progress reporting for second-pass integration for weak peaks --- .../plugins/algorithms/IntegratePeaksShoeboxTOF.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py index 2fcb2c8f34a4..7c1a3fc1fbcf 100644 --- a/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py +++ b/Framework/PythonInterface/plugins/algorithms/IntegratePeaksShoeboxTOF.py @@ -404,7 +404,9 @@ def PyExec(self): set_peak_intensity(peak, intens, sigma, do_lorz_cor) if len(ipks_strong): + prog_reporter.resetNumSteps(int(len(weak_peaks_list)), start=0.0, end=1.0) for weak_pk in weak_peaks_list: + prog_reporter.report("Re-integrating weak peaks") # get peak ipk = weak_pk.ipk peak = peaks.getPeak(ipk)