Beam center finder

In SANS experiments, it is essential to find the center of the scattering pattern in order to allow symmetric summation of the scattering intensity around the beam (i.e. computing a one-dimensional \(I(Q)\)). As detector panels can move, the center of the beam will not always be located at the same place on the detector panel from one experiment to the next.

Here we describe two different algorithms that can be used to determine the position of the beam center:

1. Using the center of mass of the pixel counts

2. Using an iterative refinement on a computed scattering cross-section to find the center of the scattering pattern

import scipp as sc
import scippnexus as snx
import plopp as pp
from ess import sans
from ess import isissans as isis
import ess.isissans.data  # noqa: F401
from ess.sans.types import *

Create and configure the workflow

We begin by setting creating a workflow and set some parameters, as well as the name of the data file we wish to use.

workflow = isis.sans2d.Sans2dTutorialWorkflow()
# For real data use:
# workflow = isis.sans2d.Sans2dWorkflow()
workflow[Filename[SampleRun]] = isis.data.sans2d_tutorial_sample_run()
workflow[isis.sans2d.LowCountThreshold] = sc.scalar(-1, unit="counts")

Downloading file 'SANS2D00063114.nxs.h5' from 'https://public.esss.dk/groups/scipp/ess/sans2d/1/SANS2D00063114.nxs.h5' to '/home/runner/.cache/ess/sans2d/1'.

Masking bad pixels

We create a quick image of the data (summing along the tof dimension) to inspect its contents. We see a diffuse scattering pattern, centered around a dark patch with an arm to the north-east; this is the sample holder. It is clear that the sample and the beam are not in the center of the panel, which is marked by the black cross

workflow[BeamCenter] = sc.vector([0, 0, 0], unit='m')
raw = workflow.compute(DetectorData[SampleRun])['spectrum', :61440]

p = isis.plot_flat_detector_xy(raw.hist(), norm='log')
p.ax.plot(0, 0, '+', color='k', ms=10)
p

The scattering pattern is circularly symmetric around the center of the beam. Because the beam is not in the center of the panel, different regions of the panel are covering different radial ranges.

When searching for the center of the beam, it is important to remove any such bias that would skew the computed center position. We basically would like to add a circular mask around the center, to ensure the same radial range is reached in all directions.

Unfortunately, we do not know the center, because that is what we are trying to compute. We can however use the fact that the scattering pattern is symmetric, and that the intensity is higher close to the beam while lower towards the edges of the panel.

Adding a mask, based on the pixel neutron counts (integrated along the tof dimension), will yield a circular mask around the beam center. Because the sample holder is highly absorbent to neutrons, such a mask will also mask out the sample holder, making it a very simple but effective way of masking.

Note

In general, for a fully normalized data reduction, masks based on counts are not recommended as they can remove true regions of zero counts and not just artifacts. They should only be used in controlled cases, where the end result does not require a correctly normalized intensity.

workflow[BeamCenter] = sc.vector(value=[0, 0, 0], unit='m')
masked = workflow.compute(MaskedData[SampleRun])['spectrum', :61440].copy()
masked.masks['low_counts'] = masked.hist().data < sc.scalar(80.0, unit='counts')

p = isis.plot_flat_detector_xy(masked.hist(), norm='log')
p.ax.plot(0, 0, '+', color='k', ms=10)
p

Method 1: center-of-mass calculation

The first method we will use to compute the center of the beam is to calculate the center-of-mass of the pixels, using the integrated counts along the time-of-flight dimension as the weights of the pixel positions.

We can visualize the pipeline that is used to compute the center-of-mass:

# The center-of-mass approach is based on the MaskedData
workflow.visualize(MaskedData[SampleRun])

We use the workflow to fill in the required arguments and call the function:

com = sans.beam_center_finder.beam_center_from_center_of_mass(workflow)
com

scipp.Variable (336 Bytes)

• • ()
    vector3
    m
    [ 0.09011603 -0.08263932 0. ]

    Values:
    array([ 0.09011603, -0.08263932,  0.        ])

We can now update our previous figure with the new position of the beam center (pink dot), which is clearly in the center of the beam/sample holder.

xc = com.fields.x
yc = com.fields.y
p.ax.plot(xc.value, yc.value, 'o', color='magenta', mec='lightgray', ms=6)
p

Method 2: computing \(I(Q)\) inside 4 quadrants

The procedure is the following:

1. Divide the panel into 4 quadrants

2. Compute \(I(Q)\) inside each quadrant and compute the residual difference between all 4 quadrants

3. Iteratively move the center position and repeat 1. and 2. until all 4 \(I(Q)\) curves lie on top of each other

For this, we need to set-up a fully-fledged workflow that can compute \(I(Q)\), which requires some additional parameters (see the SANS2D reduction workflow for more details).

Note

In the full \(I(Q)\) reduction, there is a term \(D(\lambda)\) in the normalization called the “direct beam” which gives the efficiency of the detectors as a function of wavelength. Because finding the beam center is required to compute the direct beam in the first place, we do not include this term in the computation of \(I(Q)\) for finding the beam center. This changes the shape of the \(I(Q)\) curve, but since it changes it in the same manner for all \(\phi\) angles, this does not affect the results for finding the beam center.

workflow = isis.sans2d.Sans2dTutorialWorkflow()
# For real data use:
# workflow = isis.sans2d.Sans2dWorkflow()
workflow.insert(isis.io.transmission_from_sample_run)

We set the missing input parameters:

workflow[WavelengthBins] = sc.linspace(
    'wavelength', start=2.0, stop=16.0, num=141, unit='angstrom'
)
workflow[Filename[EmptyBeamRun]] = isis.data.sans2d_tutorial_empty_beam_run()

workflow[NeXusMonitorName[Incident]] = 'monitor2'
workflow[NeXusMonitorName[Transmission]] = 'monitor4'

workflow[isis.SampleOffset] = sc.vector([0.0, 0.0, 0.053], unit='m')
workflow[isis.MonitorOffset[Transmission]] = sc.vector([0.0, 0.0, -6.719], unit='m')

workflow[CorrectForGravity] = True
workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound
workflow[DirectBeam] = None
workflow[isis.sans2d.LowCountThreshold] = sc.scalar(-1, unit="counts")

Downloading file 'SANS2D00063091.nxs.h5' from 'https://public.esss.dk/groups/scipp/ess/sans2d/1/SANS2D00063091.nxs.h5' to '/home/runner/.cache/ess/sans2d/1'.

Finally, we set the data to be used, including overriding with the new mask defined earlier:

workflow[Filename[SampleRun]] = isis.data.sans2d_tutorial_sample_run()
workflow[BeamCenter] = sc.vector([0, 0, 0], unit='m')
detector = workflow.compute(NeXusComponent[snx.NXdetector, SampleRun]).copy()
detector['data'] = detector['data']['spectrum', :61440].assign_masks(masked.masks)
workflow[NeXusComponent[snx.NXdetector, SampleRun]] = detector

The division of the panel pixels into 4 quadrants, as well as the iterative procedure to maximize the overlap between the computed intensities, is all performed internally by the beam_center_from_iofq provider (see further details below).

We can thus compute the beam center in the same way as before:

q_bins = sc.linspace('Q', 0.02, 0.25, 71, unit='1/angstrom')
workflow[BeamCenter] = sc.vector([0, 0, 0], unit='m')
iofq_center = sans.beam_center_finder.beam_center_from_iofq(
    workflow=workflow, q_bins=q_bins
)
iofq_center

scipp.Variable (336 Bytes)

• • ()
    vector3
    m
    [ 0.09266183 -0.08203386 0. ]

    Values:
    array([ 0.09266183, -0.08203386,  0.        ])

We, once again, show the location of the computed center (pink dot) on the detector panel:

p = isis.plot_flat_detector_xy(masked.hist(), norm='log')
p.ax.plot(0, 0, '+', color='k', ms=10)
p.ax.plot(
    iofq_center.value[0],
    iofq_center.value[1],
    'o',
    color='magenta',
    mec='lightgray',
    ms=6,
)
p

Finally, we can compare the values from the two methods, which should be almost identical:

print('Center-of-mass:', com.value, '\nI(Q):          ', iofq_center.value)

Center-of-mass: [ 0.09011603 -0.08263932  0.        ]
I(Q):           [ 0.09266183 -0.08203386  0.        ]

Detailed description of method 2

In the remainder of this notebook, we will describe in more detail what is done internally for method 2 in the essans module.

The user does not need to understand all the details of the implementation, the information is kept here for completeness.

Step 1: divide the panel into 4 quadrants

We divide the panel into 4 quadrants. Panels a very commonly rectangular, and the best way to ensure that each quadrant has approximately the same number of pixels is to make a vertical and a horizontal cut:

p = isis.plot_flat_detector_xy(masked.hist(), norm='log')
p.ax.axvline(0, color='cyan')
p.ax.axhline(0, color='cyan')
p.ax.plot(0, 0, '+', color='k', ms=10)
dx = 0.25
style = dict(ha='center', va='center', color='w')  # noqa: C408
p.ax.text(dx, dx, 'North-East', **style)
p.ax.text(-dx, dx, 'North-West', **style)
p.ax.text(dx, -dx, 'South-East', **style)
p.ax.text(-dx, -dx, 'South-West', **style)
p

Step 2: compute \(I(Q)\) inside each quadrant

We define several quantities which are required to compute \(I(Q)\):

workflow = workflow.copy()
workflow[QBins] = q_bins
workflow[ReturnEvents] = False
workflow[DimsToKeep] = ()
workflow[WavelengthMask] = None
workflow[WavelengthBands] = None
kwargs = dict(  # noqa: C408
    workflow=workflow,
    detector=detector['data'],
    norm=workflow.compute(CleanDirectBeam),
)

We now use a function internal to the esssans module compute \(I(Q)\) inside each quadrant:

from ess.sans.beam_center_finder import _iofq_in_quadrants

results = _iofq_in_quadrants(
    xy=[0, 0],
    **kwargs,
)

We can plot on the same figure all 4 \(I(Q)\) curves for each quadrant.

pp.plot(results, norm='log')

As we can see, the overlap between the curves from the 4 quadrants is mediocre. We will now use an iterative procedure to improve our initial guess, until a good overlap between the curves is found.

Step 3: iteratively maximize the overlap between the I(Q) curves

We first define a cost function, which gives us an idea of how good the overlap is:

\[\text{cost} = \frac{\sum_{Q}\sum_{i=1}^{i=4} \overline{I}(Q)\left(I(Q)_{i} - \overline{I}(Q)\right)^2}{\sum_{Q}\overline{I}(Q)} ~,\]

where \(\overline{I}(Q)\) is the mean intensity of the 4 quadrants (represented by \(i\)) as a function of \(Q\). This is basically a weighted mean of the square of the differences between the \(I(Q)\) curves in the 4 quadrants with respect to \(\overline{I}(Q)\), and where the weights are \(\overline{I}(Q)\).

Next, we iteratively minimize the computed cost (this is using Scipy’s optimize.minimize utility internally; see here for more details).

from scipy.optimize import minimize

# The minimizer works best if given bounds, which are the bounds of our detector panel
x = masked.coords['position'].fields.x
y = masked.coords['position'].fields.y
bounds = [(x.min().value, x.max().value), (y.min().value, y.max().value)]

res = minimize(
    sans.beam_center_finder._cost,
    x0=[0, 0],
    args=tuple(kwargs.values()),
    bounds=bounds,
    method='Powell',
    tol=0.01,
)

res

/home/runner/work/esssans/esssans/.tox/docs/lib/python3.10/site-packages/scipy/optimize/_optimize.py:2290: RuntimeWarning: invalid value encountered in scalar subtract
  p = (xf - fulc) * q - (xf - nfc) * r
/home/runner/work/esssans/esssans/.tox/docs/lib/python3.10/site-packages/scipy/optimize/_optimize.py:2291: RuntimeWarning: invalid value encountered in scalar subtract
  q = 2.0 * (q - r)

 message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 0.3043606
       x: [ 9.161e-02 -8.152e-02]
     nit: 3
   direc: [[ 1.000e+00  0.000e+00]
           [ 0.000e+00  1.000e+00]]
    nfev: 57

Once the iterations completed, the returned object contains the best estimate for the beam center:

res.x

array([ 0.09161094, -0.08152107])

We can now feed this value again into our iofq_in_quadrants function, to inspect the \(Q\) intensity in all 4 quadrants:

results = _iofq_in_quadrants(
    xy=[res.x[0], res.x[1]],
    **kwargs,
)

pp.plot(results, norm='log')

The overlap between the curves is excellent, allowing us to safely perform an azimuthal summation of the counts around the beam center.

Note

The result obtained just above is slightly different from the one obtained earlier using the workflow.

This is because in our example, we used x=0, y=0 as our initial guess, while the workflow uses an additional optimization where it first computes a better initial guess using method 1 (center-of-mass). This allows it to converge faster, with fewer iterations, and produce a slightly more accurate result.