Model fitting to small-angle neutron scattering experimental data for archaeology

Context

This project is part of an archaeological effort to identify ancient pottery-forming techniques (coil building, wheel shaping, wheel throwing) non-destructively and quantitatively, by studying the pores and grains embedded in the clay. A sherd is probed with a beam of cold neutrons, which scatter off these inclusions. The resulting 2D scattering pattern, recorded by a detector, carries information about the average shape and orientation of these inclusions, and therefore about the forming technique. The general-purpose fitting software previously used requires a user-supplied starting guess for every sample and can take minutes to converge, which does not scale to a project with hundreds of measurements.

The goal is to develop an automated procedure for analysing 2D experimental data from small-angle neutron scattering measurements of these samples. We specifically focus on the anisotropic inverse radial power model. We develop a custom preprocessing procedure followed by a maximum-likelihood fitting one to refine the parameters and compute uncertainties. This allows to reduce analysis times from several minutes to seconds per data file.

2022 internship at Centre for Energy Research in Budapest

Under a rough-surface, single-scattering hypothesis, an isotropic reference material's scattering intensity follows a similarly isotropic power-law decay. Through a Fourier transform, we show that the scattering figure of a real material can be approximated as a rotated and anisotropically scaled version of that reference. This model of the 2D scattering figure is parameterised by an intensity, roughness exponent, background, tilting angle, isotropy, and figure centre. Detector counts are modelled as Gaussian around this intensity (large-count limit of the underlying Bernoulli scattering events per pixel). The parameters are estimated by minimising the resulting χ² distance (maximum-likelihood fit).

Because this χ² function is non-convex, we design a fast preprocessing step. A weighted principal component analysis restricted to the high-count pixels gives an initial estimate of the centre, tilt and isotropy. A log-log linear regression then gives the intensity and roughness exponent. We use this as a starting point for a standard local optimiser (SQP, via the NLopt library). This methodology is not robust to measurement artefacts or high background noise, but converges correctly on some reference datasets. This reduces per-file processing time from minutes (with the previous software) to seconds.

The whole procedure is implemented in an initial version of the software, using C++14 with Qt 5.15.

Continued development

The methodology is being refined to improve robustness against experimental noise and to better fit the chosen model. The current version of the software is developed in Rust using iced.