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 the beam center finding algorithm, which uses a combination of a center-of-mass calculation and an iterative refinement on a computed scattering cross-section to find the center of the scattering pattern.

[1]:

import scipp as sc
from scipp.constants import g
from ess import loki, sans
from ess.sans import beam_center_finder as bcf
from ess.logging import configure_workflow
import plopp as pp

[2]:

logger = configure_workflow('sans_beam_center_finder', filename='sans.log')

We begin by setting some parameters relevant to the current sample.

[3]:

# Include effects of gravity?
gravity = True

# Wavelength binning
wavelength_bins = sc.linspace(
    dim='wavelength', start=2.0, stop=16.0, num=141, unit='angstrom'
)

# Define Q binning
q_bins = sc.linspace('Q', 0.02, 0.3, 71, unit='1/angstrom')

# Define coordinate transformation graph
graph = sans.conversions.sans_elastic(gravity=gravity)

Next we load the data files for the sample and direct runs:

[4]:

sample = loki.io.load_sans2d(filename=loki.data.get_path('SANS2D00063114.hdf5'))
direct = loki.io.load_sans2d(filename=loki.data.get_path('SANS2D00063091.hdf5'))
# Add X, Y coordinates
sample.coords['x'] = sample.coords['position'].fields.x
sample.coords['y'] = sample.coords['position'].fields.y
# Add gravity coordinate
sample.coords["gravity"] = sc.vector(value=[0, -1, 0]) * g

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 red dot

[5]:

image = sample.sum('tof').copy().hist(y=120, x=128)
p = image.plot(norm='log', aspect='equal')
p.ax.plot(0, 0, 'o', color='red', ms=5)
p

[5]:

../../_images/techniques_sans_sans-beam-center-finder_8_0.svg

To avoid skew in future comparisons of integrated intensities between the different regions of the detector panel, we mask out the sample holder, using a low-counts threshold. This also masks out the edges of the panel, which show visible artifacts. We also mask out a region in the bottom right corner where a group of hot pixels is apparent. Finally, there is a single hot pixel in the detector on the right edge of the panel with counts in excess of 1000, which we also remove.

[6]:

summed = sample.sum('tof')

low_counts = summed.data < sc.scalar(70, unit='counts')
high_counts = summed.data > sc.scalar(1000, unit='counts')
lower_right = (sample.coords['x'] > sc.scalar(0.35, unit='m')) & (
    sample.coords['y'] < sc.scalar(-0.4, unit='m')
)

sample.masks['low_counts'] = low_counts
sample.masks['high_counts'] = high_counts
sample.masks['lower_right'] = lower_right

We look at the image again, to verify we have masked the desired regions.

[7]:

image = sample.sum('tof').copy().hist(y=120, x=128)
p = image.plot(norm='log', aspect='equal')
p.ax.plot(0, 0, 'o', color='red', ms=5)
p

[7]:

../../_images/techniques_sans_sans-beam-center-finder_12_0.svg

Description of the procedure#

The procedure to determine the precise location of the beam center is the following:

obtain an initial guess by computing the center-of-mass of the pixels, weighted by the counts on each pixel
from that initial guess, divide the panel into 4 quadrants
compute \(I(Q)\) inside each quadrant and compute the residual difference between all 4 quadrants
iteratively move the centre position and repeat 2. and 3. until all 4 \(I(Q)\) curves lie on top of each other

Initial guess: center-of-mass calculation#

Computing the center-of-mass is straightforward using the vector position coordinate.

[8]:

com = bcf.center_of_mass(sample)

# We compute the shift between the incident beam direction and the center-of-mass
incident_beam = sample.transform_coords('incident_beam', graph=graph).coords[
    'incident_beam'
]
n_beam = incident_beam / sc.norm(incident_beam)
com_shift = com - sc.dot(com, n_beam) * n_beam
com_shift

[8]:

scipp.Variable (336 Bytes)

()

vector3

[ 0.0642058 -0.05088138 0. ]

Values:
array([ 0.0642058 , -0.05088138,  0.        ])

We can now plot the center-of-mass on the same image as before:

[9]:

xc = com_shift.fields.x
yc = com_shift.fields.y
p = image.plot(norm='log', aspect='equal')
p.ax.plot(xc.value, yc.value, 'o', color='red', ms=5)
p

[9]:

../../_images/techniques_sans_sans-beam-center-finder_16_0.svg

Making 4 quadrants#

We divide the panel into 4 quadrants.

[10]:

p = image.plot(norm='log', aspect='equal')
p.ax.plot(xc.value, yc.value, 'o', color='red', ms=5)
p.ax.axvline(xc.value, color='cyan')
p.ax.axhline(yc.value, color='cyan')
dx = 0.25
p.ax.text(xc.value + dx, yc.value + dx, 'North-East', ha='center', va='center')
p.ax.text(xc.value - dx, yc.value + dx, 'North-West', ha='center', va='center')
p.ax.text(xc.value - dx, yc.value - dx, 'South-East', ha='center', va='center')
p.ax.text(xc.value + dx, yc.value - dx, 'South-West', ha='center', va='center')
p

[10]:

../../_images/techniques_sans_sans-beam-center-finder_18_0.svg

Converting to \(Q\) inside each quadrant#

We now perform a full:math:`^1 <#footnote1>`__ \(I(Q)\) reduction (see here for more details) inside each quadrant. The reduction involves computing a normalizing term which, for the most part, does not depend on pixel positions. We therefore compute this once, before starting iterations to refine the position of the center.

Compute \(I(Q)\) inside the 4 quadrants#

We begin by defining several parameters which are required to compute \(I(Q)\).

[13]:

# Define 4 phi bins
pi = sc.constants.pi.value
phi_bins = sc.linspace('phi', -pi, pi, 5, unit='rad')

# Name the quadrants
quadrants = ['south-west', 'south-east', 'north-east', 'north-west']

# Define a wavelength range to use
wavelength_range = sc.concat(
    [wavelength_bins.min(), wavelength_bins.max()], dim='wavelength'
)

We now use a function which will apply the center offset to the pixel coordinates, and compute \(I(Q)\) inside each quadrant.

[14]:

grouped = sans.beam_center_finder.iofq_in_quadrants(
    xy=[com_shift.fields.x.value, com_shift.fields.y.value],
    sample=sample,
    norm=norm,
    graph=graph,
    q_bins=q_bins,
    gravity=gravity,
    wavelength_range=wavelength_range,
)

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

[15]:

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

[15]:

../../_images/techniques_sans_sans-beam-center-finder_28_0.svg

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

For this, we first define a cost function sans.beam_center_finder.cost, 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).

[16]:

# The minimizer works best if given bounds, which are the bounds of our detector panel
x = sample.coords['position'].fields.x
y = sample.coords['position'].fields.y

res = bcf.minimize(
    sans.beam_center_finder.cost,
    x0=[com_shift.fields.x.value, com_shift.fields.y.value],
    args=(sample, norm, graph, q_bins, gravity, wavelength_range),
    bounds=[(x.min().value, x.max().value), (y.min().value, y.max().value)],
)
res

[16]:

       message: Optimization terminated successfully.
       success: True
        status: 0
           fun: 0.9083893895149231
             x: [ 9.330e-02 -8.079e-02]
           nit: 20
          nfev: 37
 final_simplex: (array([[ 9.330e-02, -8.079e-02],
                       [ 9.347e-02, -8.129e-02],
                       [ 9.242e-02, -8.121e-02]]), array([ 9.084e-01,  9.100e-01,  9.882e-01]))

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

[17]:

res.x

[17]:

array([ 0.09329701, -0.08079138])

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

[18]:

grouped = sans.beam_center_finder.iofq_in_quadrants(
    xy=[res.x[0], res.x[1]],
    sample=sample,
    norm=norm,
    graph=graph,
    q_bins=q_bins,
    gravity=gravity,
    wavelength_range=wavelength_range,
)

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

[18]:

../../_images/techniques_sans_sans-beam-center-finder_34_0.svg

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

As a consistency check, we plot the refined beam center position onto the detector panel image:

[19]:

p = sample.sum('tof').copy().hist(y=120, x=128).plot(norm='log', aspect='equal')
p.ax.plot(res.x[0], res.x[1], 'o', color='red', ms=5)
p

[19]:

../../_images/techniques_sans_sans-beam-center-finder_36_0.svg

Footnotes#

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.

Beam center finder

Contents