BornAgain ≻ 24 ≻ Documentation ≻ Reference ≻ Fitting ≻ Examples ≻ Experiment at GALAXI

Experiment at GALAXI

This is an example of a real data fit. We use our own measurements performed at the laboratory diffractometer GALAXI in Forschungszentrum Jülich.

Real-space model

Historical fit (BornAgain 1.19)

  • The sample represents a 4 layer system (substrate, teflon, hmdso and air) with Ag nanoparticles placed inside the hmdso layer on top of the teflon layer.
  • The nanoparticles have a broad log-normal size distribution, modeled by a Mixture of spheres with different radii.
  • All three interfaces have self-affine roughness. A CommonDepthCrosscorrelation vertically correlates the upper interfaces with the interface below, producing the roughness-interference contribution alongside the particle scattering.
  • The detector is created to represent the PILATUS detector from the experiment, with angular ranges anchored to the direct-beam position.
  • The experimental data and detector are cropped to the fit window. A detector bitmap mask then excludes the specular beam; masked pixels appear as NaN in the simulation result and are not computed.
  • The real data is loaded from a tiff file into a numpy array.
  • The residual function compares the data with the simulation under a Poisson-like weighting. The mean particle radius and distance retain their reported values. The fit varies the width of the particle-size distribution, two interface RMS roughnesses, an effective intensity normalization, and a constant background using least-squares minimization.

Experiment description

To successfully simulate and fit results of some real experiment it is important to have

  • A good guess about the sample structure and the initial values of the sample parameters.
  • Full information about the instrument geometry: size and angular range of the detector.
  • A 2D numpy array containing the intensities measured in the detector channels.

Experiment

As an example we will use our own measurements performed at the laboratory diffractometer GALAXI in Forschungszentrum Jülich.

Our sample represents a 4-layer system (substrate, teflon, hmdso and air) with Ag nanoparticles embedded in the hmdso layer, on top of the teflon layer. The PILATUS 1M detector was placed at a distance of 1730 mm from the sample.

The results of the measurement are represented by the intensity image taken in certain conditions (beam wavelength, grazing angle, detector angular range) and stored in a 32-bit tiff file. To be able to fit these data we have to

  • prepare a description of the simulation
  • load the experimental data in BornAgain’s fitting engine

Preparing the simulation description

From the experimental setup we know the pixel size, the sample-detector distance, and the point where the direct beam hits the detector. Two helper functions convert positions and pixel edges on the flat detector to scattering angles in radians:

def detector_angle(position, beam_position, detector_distance, offset=0):
    return np.arctan2(position - beam_position, detector_distance) - offset


def detector_angle_range(pixel_centers, pixel_size, beam_position,
                         detector_distance, offset=0):
    lower_edge = pixel_centers[0] - pixel_size/2
    upper_edge = pixel_centers[-1] + pixel_size/2
    return (detector_angle(lower_edge, beam_position, detector_distance,
                           offset),
            detector_angle(upper_edge, beam_position, detector_distance,
                           offset))

The masks and detector are then constructed explicitly from the calibrated coordinates:

full_nx = 981
full_ny = 1043
pixel_size = 0.172  # in mm
detector_distance = 1730  # in mm
beam_x_pixel_pos, beam_y_pixel_pos = 597.1, 323.4  # from lower left
alpha_i = 0.463*deg

x_pixels = (np.arange(full_nx) + 0.5)*pixel_size
y_pixels = (np.arange(full_ny) + 0.5)*pixel_size
x_grid, y_grid = np.meshgrid(x_pixels, y_pixels)
inside_window = ((x_grid > 85) & (x_grid < 120)
                 & (y_grid > 70) & (y_grid < 92))
window_mask = np.logical_not(inside_window)
fit_x = ba_fit.crop_by_mask(x_grid, window_mask)
fit_y = ba_fit.crop_by_mask(y_grid, window_mask)

beam_x_pos = beam_x_pixel_pos*pixel_size
beam_y_pos = beam_y_pixel_pos*pixel_size
phi_min, phi_max = detector_angle_range(
    fit_x[0, :], pixel_size, beam_x_pos, detector_distance)
alpha_min, alpha_max = detector_angle_range(
    fit_y[:, 0], pixel_size, beam_y_pos, detector_distance, alpha_i)

specular_beam_mask = ((fit_x > 101.9) & (fit_x < 103.7)
                      & (fit_y > 82.1) & (fit_y < 85.2))

cropped_ny, cropped_nx = fit_x.shape
detector = ba.SphericalDetector(
    cropped_nx, phi_min, phi_max,
    cropped_ny, alpha_min, alpha_max)
detector.setMask(specular_beam_mask)

The fit window and specular-beam exclusion remain in millimeters, matching their experimental specification. ba_fit.crop_by_mask follows the detector mask convention: True marks excluded pixels. It returns the part of the data inside the smallest rectangle containing every unmasked pixel.

The experimental array is cropped to the fit window before fitting. A spherical detector with the same number of pixels and angular boundaries then approximates this part of the flat detector to a fraction of a pixel.

Masking detector pixels

After cropping, a boolean bitmap mask excludes the specular beam. The mask shape is (n_alpha, n_phi) and True marks excluded pixels.

The final simulation setup looks as follows:

beam = ba.Beam(P["intensity"], 1.34*angstrom, alpha_i)
sample = get_sample(P)
simulation = ba.ScatteringSimulation(beam, sample, detector)

During the fit, masked detector pixels are stored as NaN in the simulation result. The residual function therefore only reads the unmasked pixels, and additionally skips the negative dead-pixel markers of the PILATUS image.

Importing the real data using Fabio library

The Fabio library provides a convenient way to import experimental data in the form of a numpy array.

import fabio

img = fabio.open("galaxi_data.tif.gz")
data = np.flipud(img.data.astype("float64"))
data = ba_fit.crop_by_mask(data, window_mask)

The full image has shape (1043, 981). After np.flipud, the fit window is selected by the boolean mask and has the same shape as the cropped detector, (128, 204). Note that fabio returns the image with row 0 at the top detector row, whereas BornAgain expects row 0 at the smallest scattering angle.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

#!/usr/bin/env python3
# /// script
# requires-python = ">=3.10"
# dependencies = ["bornagain>=25,<26", "fabio", "lmfit"]
# ///
"""
Fitting experimental data: spherical nanoparticles with size distribution
in a four-layer system (experiment at GALAXI).

Demonstrates how to load a detector image with fabio, anchor the
detector geometry to the direct-beam position, crop the data and
detector to the fit window, and compare a polydisperse sample model
to the data by fitting the particle-size distribution width, interface
roughnesses, intensity normalization, and constant background.
"""
import os
import bornagain as ba
ba.require_versions("bornagain>=25,<26")
from bornagain import angstrom, ba_fit, ba_fitmonitor, ba_io, deg, nm, R3
from matplotlib import pyplot as plt
import lmfit
import numpy as np

datadir = ba_io.data_dir()


def get_mixture(ag_mat, radius, sigma, hmdso_thickness):
    """
    Creates a log-normal mixture of Ag spheres that fit in the HMDSO layer.
    """
    nsizes = 20
    r_min = radius*np.exp(-2*sigma)
    r_max = min(radius*np.exp(2*sigma), hmdso_thickness/2)
    radii = np.linspace(r_min, r_max, nsizes)
    weights = np.exp(-np.log(radii/radius)**2/(2*sigma**2))/radii
    mixture = ba.Mixture()  # normalizes the weights
    for r, w in zip(radii, weights):
        sphere = ba.Particle(ag_mat, ba.Sphere(r))
        sphere.translate(R3(0, 0, -hmdso_thickness))
        mixture.addParticle(sphere, w)
    return mixture


def get_sample(P):
    """
    Creates the layered GALAXI sample for the fit parameters P.
    """
    radius = P["radius"]*nm
    sigma = P["sigma"]
    distance = P["distance"]*nm
    disorder = 10.5*nm
    kappa = 17.5
    ptfe_thickness = 22.1*nm
    hmdso_thickness = 18.5*nm

    # defining materials
    vacuum = ba.Vacuum()
    si_color = (0.30, 0.62, 0.86)
    si_mat = ba.RefractiveMaterial("Si", si_color, 5.7816e-6, 1.0229e-7)
    ag_color = (0.86, 0.24, 0.18)
    ag_mat = ba.RefractiveMaterial("Ag", ag_color, 2.2475e-5, 1.6152e-6)
    ptfe_color = (0.93, 0.72, 0.25)
    ptfe_mat = ba.RefractiveMaterial("PTFE", ptfe_color, 5.20509e-6, 1.9694e-8)
    hmdso_color = (0.25, 0.65, 0.35)
    hmdso_mat = ba.RefractiveMaterial("HMDSO", hmdso_color, 2.0888e-6, 1.3261e-8)

    # collection of particles with size distribution
    mixture = get_mixture(ag_mat, radius, sigma, hmdso_thickness)

    # interference function
    layout = ba.RadialParacrystal(mixture, distance, 1e6*nm)
    layout.setKappa(kappa)
    layout.setDomainSize(2e4*nm)
    profile = ba.Profile1DGauss(disorder)
    layout.setProbabilityDistribution(profile)

    vertical_correlation_depth = 100*nm
    hmdso_roughness_rms = P["hmdso_rms"]*nm
    ptfe_roughness_rms = P["ptfe_rms"]*nm
    # The unmeasured PTFE/Si interface uses the PTFE roughness scale.
    substrate_roughness_rms = ptfe_roughness_rms
    roughness_hurst = 0.3
    lateral_correlation_length = 5*nm

    # Self-affine roughness on all interfaces. The upper two interfaces
    # are vertically correlated with the interface below.
    hmdso_autocorr = ba.SelfAffineFractalModel(
        hmdso_roughness_rms, roughness_hurst,
        lateral_correlation_length)
    ptfe_autocorr = ba.SelfAffineFractalModel(
        ptfe_roughness_rms, roughness_hurst,
        lateral_correlation_length)
    substrate_autocorr = ba.SelfAffineFractalModel(
        substrate_roughness_rms, roughness_hurst,
        lateral_correlation_length)
    transient = ba.TanhTransient()
    crosscorr = ba.CommonDepthCrosscorrelation(vertical_correlation_depth)
    hmdso_roughness = ba.Roughness(
        hmdso_autocorr, transient, crosscorr)
    ptfe_roughness = ba.Roughness(
        ptfe_autocorr, transient, crosscorr)
    substrate_roughness = ba.Roughness(substrate_autocorr, transient)

    # layers
    vacuum_layer = ba.Layer(vacuum)
    hmdso_layer = ba.Layer(hmdso_mat, hmdso_thickness, hmdso_roughness)
    hmdso_layer.deposit2D(layout)
    ptfe_layer = ba.Layer(ptfe_mat, ptfe_thickness, ptfe_roughness)
    substrate_layer = ba.Layer(si_mat, substrate_roughness)

    # assembling sample
    sample = ba.Sample()
    sample.addLayer(vacuum_layer)
    sample.addLayer(hmdso_layer)
    sample.addLayer(ptfe_layer)
    sample.addLayer(substrate_layer)

    return sample


def detector_angle(position, beam_position, detector_distance, offset=0):
    """
    Converts a position on a flat detector to a scattering angle in radians.
    """
    return np.arctan2(position - beam_position, detector_distance) - offset


def detector_angle_range(pixel_centers, pixel_size, beam_position,
                         detector_distance, offset=0):
    """
    Returns the angular range covered by the given flat-detector pixels.
    """
    lower_edge = pixel_centers[0] - pixel_size/2
    upper_edge = pixel_centers[-1] + pixel_size/2
    return (detector_angle(lower_edge, beam_position, detector_distance,
                           offset),
            detector_angle(upper_edge, beam_position, detector_distance,
                           offset))


def get_simulation(P, detector, alpha_i):
    """
    Creates a GISAS simulation for the fit parameters P.
    """
    wavelength = 1.34*angstrom
    beam = ba.Beam(P["intensity"], wavelength, alpha_i)
    sample = get_sample(P)
    return ba.ScatteringSimulation(beam, sample, detector)


def load_data(filename, window_mask):
    """
    Loads, orients, and crops the experimental detector image.

    Fabio returns row 0 at the top; BornAgain expects row 0 at the
    smallest scattering angle.
    """
    import fabio
    filepath = os.path.join(datadir, filename)
    data = np.flipud(fabio.open(filepath).data.astype(float))
    return ba_fit.crop_by_mask(data, window_mask)


def datafield_from_array(frame, array):
    """
    Wraps a NumPy array in a Datafield for plotting.
    """
    return ba.Datafield(frame, np.asarray(array).ravel().tolist())


if __name__ == '__main__':
    # Detector setup as given by the instrument responsible.
    full_nx = 981
    full_ny = 1043
    pixel_size = 0.172  # in mm
    detector_distance = 1730  # in mm
    beam_x_pixel_pos, beam_y_pixel_pos = 597.1, 323.4  # from lower left
    alpha_i = 0.463*deg

    # Pixel-center coordinates and fit window, in mm on the detector.
    x_pixels = (np.arange(full_nx) + 0.5)*pixel_size
    y_pixels = (np.arange(full_ny) + 0.5)*pixel_size
    x_grid, y_grid = np.meshgrid(x_pixels, y_pixels)
    inside_window = ((x_grid > 85) & (x_grid < 120)
                     & (y_grid > 70) & (y_grid < 92))
    window_mask = np.logical_not(inside_window)
    fit_x = ba_fit.crop_by_mask(x_grid, window_mask)
    fit_y = ba_fit.crop_by_mask(y_grid, window_mask)

    # Convert the selected pixel edges to the spherical-detector angles.
    beam_x_pos = beam_x_pixel_pos*pixel_size
    beam_y_pos = beam_y_pixel_pos*pixel_size
    phi_min, phi_max = detector_angle_range(
        fit_x[0, :], pixel_size, beam_x_pos, detector_distance)
    alpha_min, alpha_max = detector_angle_range(
        fit_y[:, 0], pixel_size, beam_y_pos, detector_distance, alpha_i)

    # Exclude the specular beam from the cropped detector.
    specular_beam_mask = ((fit_x > 101.9) & (fit_x < 103.7)
                          & (fit_y > 82.1) & (fit_y < 85.2))

    cropped_ny, cropped_nx = fit_x.shape
    detector = ba.SphericalDetector(
        cropped_nx, phi_min, phi_max,
        cropped_ny, alpha_min, alpha_max)
    detector.setMask(specular_beam_mask)

    data = load_data("scatter2d/galaxi_data.tif.gz", window_mask)
    # usable pixels: outside the specular beam, no detector gaps (-1)
    usable = (np.logical_not(specular_beam_mask)
              & np.isfinite(data) & (data >= 0))
    flat_exp_values = data[usable]
    sim_result = None  # latest simulation, shared with the plot callback

    def residuals(P):
        """
        Returns Poisson-weighted residuals for usable pixels.
        """
        global sim_result
        p = P.valuesdict()
        sim_result = get_simulation(p, detector, alpha_i).simulate()
        flat_sim_values = sim_result.intensities()[usable]
        model = flat_sim_values + p["background"]
        return (flat_exp_values - model)/np.sqrt(np.maximum(1., model))

    observer = ba_fitmonitor.PlotterGISAS(zmin=5, zmax=1e3)
    best_result = None
    best_params = None
    best_chi2 = np.inf

    def plot_iteration(P, iteration, resid):
        """
        Updates the fit-progress plot with the best result found so far.
        """
        global best_result, best_params, best_chi2
        chi2 = float(np.sum(resid*resid))
        if sim_result is not None and chi2 < best_chi2:
            best_result = sim_result
            best_params = P.copy()
            best_chi2 = chi2
        if iteration % 50 == 0 and best_result is not None:
            display_data = np.where(usable, data, np.nan)
            exp_data = datafield_from_array(best_result.frame(),
                                            display_data)
            model_data = datafield_from_array(
                best_result.frame(),
                best_result.intensities()
                + best_params["background"].value)
            observer.plot(exp_data, model_data, best_params, best_chi2)

    P = lmfit.Parameters()
    # The mean particle radius and distance retain their reported values.
    P.add("radius", value=5.75, vary=False)  # (nm)
    P.add("sigma", value=0.4, min=0., max=3., vary=True)  # (dimensionless)
    P.add("distance", value=53.6, vary=False)  # (nm)
    # RMS roughnesses start from the values reported for the sample.
    P.add("hmdso_rms", value=1.1, min=0., max=5.)  # (nm)
    P.add("ptfe_rms", value=2.3, min=0., max=5.)  # (nm)
    # Effective normalization combines source flux and particle density.
    P.add("intensity", value=4e10, min=1e9, max=1e12)  # (a.u.)
    P.add("background", value=1., min=0., max=50.)  # (counts/pixel)

    result = lmfit.minimize(
        residuals, P, method="leastsq", iter_cb=plot_iteration)
    final_residual = residuals(result.params)
    plot_iteration(result.params, 0, final_residual)
    print(lmfit.fit_report(result))

    plt.show()
auto/Examples/fit/scatter2d/expfit_galaxi.py

Data to be fitted: galaxi_data.tif.gz