INTER psi scan

[Tom-Willemsen]

(sub-ticket of #210)

For a psi (rocking) scan on INTER, we should be able to use detector-mapping to directly calculate psi, from a single data collection, rather than scanning. This should save quite a bit of time during sample alignment relative to scanning, or manually tweaking psi.

For a 2-d detector map, we will get a very elongated 2-dimensional ellipsoid. We want to extract a 'rotation' parameter from that ellipsoid; that rotation parameter is the psi value that we need.

---

Or, in other words:

• Count once
• Get this image on the detector

• Find angle of ellipse relative to horizontal, that is equal to psi
• Change psi to make it zero
• Profit

---

Workflow

• Read in full DAE data via detector-mapping interface
• Convert pixels to x/y coordinates
• Pass to lmfit
  • Documentation about how to handle 2-d data in lmfit, including elongated 2-d ellipsoids that look very similar to what we need: https://lmfit.github.io/lmfit-py/examples/example_two_dimensional_peak.html
• Return fitted angle of the ellipsoid

Acceptance criteria

• Psi scan using the above workflow can be done on INTER
• Deployed/tested on INTER

[Tom-Willemsen]

Scientist says:

beam on Inter is not pollinated in the horizontal direction. For small angles it may work well enough. I'll think about it a bit more

Hold off implementing this until we have confirmation from scientists this is the approach they want to use

[Tom-Willemsen]

[Tom-Willemsen]

[Tom-Willemsen]

Script to fit above data

import h5py
import numpy as np
import numpy.typing as npt
import matplotlib.pyplot as plt
import scipp as sc
import lmfit
from lmfit.lineshapes import gaussian
import matplotlib


matplotlib.rcParams['figure.autolayout'] = True
matplotlib.rcParams['font.size'] = 8


PIXEL_SIZE = 0.606
DETECTOR_BANKS = 8
DETECTOR_PIXELS = 256
DETECTOR_OFFSET = 4


def model_func(x, y, amplitude, centerx, centery, width,
                 rotation):
    """
    Model for a 2-D Gaussian line profile, with rotation.

    Args:
        x: x coordinates
        y: y coordinates
        amplitude: amplitude of the Gaussian
        centerx: x value of a point P through which the line passes
        centery: y value of a point P through which the line passes
        width: width of Gaussian (sigma)
        rotation: rotation of the Gaussian (radians, clockwise positive)
    """
    # See https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_point_and_angle
    dist_from_line = np.abs(np.cos(rotation)*(centery-y) - np.sin(rotation)*(centerx - x))

    return gaussian(dist_from_line, amplitude=amplitude, sigma=width)

def load_inter_data():
    """
    Load testing data.

    In operation, this would be replaced by g.get_specdata().
    """
    with h5py.File("c:\\instrument\\tmp\\INTER00086142.nxs", "r") as h5:
        return h5["raw_data_1"]["detector_1"]["counts"][0,:,:]


def reshape_data(inp: npt.NDArray[np.float32]) -> npt.NDArray[np.float64]:
    """
    Args:
        inp (npt.NDArray[np.float64]): a 2-dimensional numpy array (detectors, time channels)

    Returns:
        A 2-dimensional numpy array (detector x, y image), summed in ToF.

    """
    data = np.astype(inp, "float64")

    data = data.sum(axis=1)  # sum along the time-of-flight dimension

    # remove monitor spectra
    data = data[DETECTOR_OFFSET:DETECTOR_OFFSET + DETECTOR_PIXELS * DETECTOR_BANKS]

    data = data.reshape((DETECTOR_BANKS, DETECTOR_PIXELS))
    data = data[::-1, ::-1]
    return data.transpose()


def rebin_and_normalize(data: npt.NDArray[np.float64]) -> sc.DataArray:
    """
    Rebin and normalize the data.

    Args:
        data (npt.NDArray[np.float64]): a 2-dimensional numpy array of (x, y) pixel
            intensities, summed in ToF.

    Returns:
        DataArray, rebinned to square pixels and normalized by the sum of intensity
            in each detector bank.

    """
    roi_start = 150 * PIXEL_SIZE
    roi_end = 200 * PIXEL_SIZE

    var = sc.DataArray(
        data=sc.array(dims=["pixel", "bank"], values=data, unit=sc.units.counts),
        coords={
            "bank": sc.array(dims=["bank"], values=[-210, -150, -90, -30, 0, 30, 90, 150, 210], unit=sc.units.mm,
                             dtype="float64"),
            "pixel": sc.linspace(dim="pixel", start=-0.303, stop=155.136 + 0.303, num=256 + 1, dtype="float64",
                                 unit=sc.units.mm),
        }
    )

    # Exclude outer banks, they don't usually count much.
    rebin_coord_bank = sc.arange(dim="bank", start=-150., stop=150., step=PIXEL_SIZE, dtype="float64", unit=sc.units.mm)
    rebin_coord_pixel = sc.arange(dim="pixel", start=roi_start, stop=roi_end, step=PIXEL_SIZE, dtype="float64",
                                  unit=sc.units.mm)

    # Rebin to square pixels
    var_rebinned = var.rebin({
        "bank": rebin_coord_bank,
        "pixel": rebin_coord_pixel,
    })

    # Normalize so all per-bank sums are equal
    return var_rebinned / var_rebinned.sum(dim="pixel")


def draw_plots(data: npt.NDArray, model: lmfit.model.Model, result: lmfit.model.ModelResult, x_cen: npt.NDArray, y_cen: npt.NDArray):
    """
    Show fit quality visually using matplotlib.
    """
    X, Y = np.meshgrid(x_cen, y_cen)

    fit = model.func(X, Y, **result.best_values)

    (_, (ax1, ax2, ax3)) = plt.subplots(3)

    ax1.imshow(data, aspect="auto", cmap="inferno", vmin=0, vmax=np.max(data))
    ax2.imshow(fit, aspect="auto", cmap="inferno", vmin=0, vmax=np.max(data))
    ax3.imshow(np.abs(fit - data), aspect="auto", cmap="inferno", vmin=0, vmax=np.max(data))

    ax1.set_title("Data")
    ax2.set_title("Fit")
    ax3.set_title("data-fit")

    plt.show()


def calculate_psi(inp: npt.NDArray[np.float32]) -> float:
    """
    Calculate psi from a detector-map.

    Parameters:
        inp: a 2-dimensional numpy array (detectors, time channels).
            This should not include the 'junk' time channel zero or spectrum zero.

    Returns:
        psi (degrees, clockwise rotation)

    """
    data = reshape_data(inp)
    da = rebin_and_normalize(data)

    data = da.values

    x = da.coords["bank"].values
    y = da.coords["pixel"].values

    x_cen = (x[:-1] + x[1:]) / 2.
    y_cen = (y[:-1] + y[1:]) / 2.

    model = lmfit.Model(model_func, independent_vars=['x', 'y'])
    params = model.make_params()
    params["rotation"].set(value=0, min=-np.pi/2, max=np.pi/2, vary=True)
    params["centerx"].set(0, vary=False)
    params["centery"].set((np.max(y) + np.min(y)) / 2, min=np.min(y), max=np.max(y), vary=True)
    params["amplitude"].set(np.max(data), min=0, vary=True)
    params["width"].set(value=PIXEL_SIZE, min=0., vary=True)

    result = model.fit(
        data,
        x=x_cen.reshape(1, len(x_cen)),
        y=y_cen.reshape(len(y_cen), 1),
        params=params,
    )

    print(result.fit_report(show_correl=False))

    draw_plots(data, model, result, x_cen, y_cen)

    # Factor of two is because image on detector is rotated
    # twice as much as the actual psi angle at the sample.
    # See diagrams from INTER scientists for further clarification
    return np.rad2deg(result.values["rotation"]) / 2.0


if __name__ == "__main__":
    d = load_inter_data()
    print(f"PSI = {calculate_psi(d):.3f} deg")

---

The above script gives a fit like this (after some normalization and rebinning steps) with a Gaussian line profile, which seemed to be fitting a bit better than a Gaussian blob:

With fitting parameters:

    amplitude:  0.56571282 +/- 0.00179698 (0.32%) (init = 0.3943265)
    centerx:    0 (fixed)
    centery:    109.279173 +/- 0.00291769 (0.00%) (init = 106.05)
    width:      0.79537996 +/- 0.00291739 (0.37%) (init = 0.606)
    rotation:  -0.01552527 +/- 3.3686e-05 (0.22%) (init = 0)

psi sould be rad2deg(rotation) / 2.0. The factor of 2 is because rotation of projection on detector is doubled relative to rotation of sample (see diagrams above).

I think this script should live in INTER's instrument scripts, it is conceptually nothing to do with bluesky and highly instrument-specific.

To test on INTER:

• Add this script to INTER's inst scripts
• Swap out the HDF file loading with a g.get_specdata call or a bluesky bps.count() plan with a PeriodSpecIntegralsReducer configured.
• Try various psi values and check that this method actually gets the correct answer for all realistic samples / values of psi (this will need input from scientists)

---

Once we have tested this, this psi alignment can be included as part of INTER's normal sample alignment/scanning scripts.