Real life fit example: Experiment at GALAXI

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

Real-space model: 
Intensity Image: 
Python Script: 

The file SampleBuilder.py contains a sample description

  • 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.
  • Sample is generated with the help of SampleBuilder, which is able to create samples depending from registered parameters
  • Nanoparticles have a broad Log-normal size distribution

 

"""
3 layers system (substrate, teflon, air).
Air layer is populated with spheres with some size distribution.
"""
import bornagain as ba
import ctypes


class MySampleBuilder(ba.IMultiLayerBuilder):
    """

    """
    def __init__(self):
        ba.IMultiLayerBuilder.__init__(self)
        self.sample = None

        # parameters describing the sample
        self.radius = ctypes.c_double(5.75*ba.nm)
        self.sigma = ctypes.c_double(0.4)
        self.distance = ctypes.c_double(53.6*ba.nm)
        self.disorder = ctypes.c_double(10.5*ba.nm)
        self.kappa = ctypes.c_double(17.5)
        self.ptfe_thickness = ctypes.c_double(22.1*ba.nm)
        self.hmdso_thickness = ctypes.c_double(18.5*ba.nm)

        # register parameters
        self.registerParameter("radius", ctypes.addressof(self.radius))
        self.registerParameter("sigma", ctypes.addressof(self.sigma))
        self.registerParameter("distance", ctypes.addressof(self.distance))
        self.registerParameter("disorder", ctypes.addressof(self.disorder))
        self.registerParameter("kappa", ctypes.addressof(self.kappa))
        self.registerParameter("tptfe", ctypes.addressof(self.ptfe_thickness))
        self.registerParameter("thmdso", ctypes.addressof(self.hmdso_thickness))

    # constructs the sample for current values of parameters
    def buildSample(self):
        # defining materials
        m_air = ba.HomogeneousMaterial("Air", 0.0, 0.0)
        m_Si = ba.HomogeneousMaterial("Si", 5.78164736e-6, 1.02294578e-7)
        m_Ag = ba.HomogeneousMaterial("Ag", 2.24749529E-5, 1.61528396E-6)
        m_PTFE = ba.HomogeneousMaterial("PTFE", 5.20508729E-6, 1.96944292E-8)
        m_HMDSO = ba.HomogeneousMaterial("HMDSO", 2.0888308E-6, 1.32605651E-8)

        # collection of particles with size distribution
        nparticles = 20
        nfwhm = 2.0
        sphere_ff = ba.FormFactorFullSphere(self.radius.value)
        # sphere_ff = ba.FormFactorTruncatedSphere(
        #    self.radius.value, self.radius.value*1.5)

        sphere = ba.Particle(m_Ag, sphere_ff)
        position = ba.kvector_t(0*ba.nm, 0*ba.nm,
                                -1.0*self.hmdso_thickness.value)
        sphere.setPosition(position)
        ln_distr = ba.DistributionLogNormal(self.radius.value, self.sigma.value)
        par_distr = ba.ParameterDistribution(
            "/Particle/FullSphere/Radius", ln_distr, nparticles, nfwhm)
        # par_distr = ba.ParameterDistribution(
        #    "/Particle/TruncatedSphere/Radius", ln_distr, nparticles, nfwhm)
        # par_distr.linkParameter("/Particle/TruncatedSphere/Height")
        part_coll = ba.ParticleDistribution(sphere, par_distr)

        # interference function
        interference = ba.InterferenceFunctionRadialParaCrystal(
            self.distance.value, 1e6*ba.nm)
        interference.setKappa(self.kappa.value)
        interference.setDomainSize(20000.0)
        pdf = ba.FTDistribution1DGauss(self.disorder.value)
        interference.setProbabilityDistribution(pdf)

        # assembling particle layout
        particle_layout = ba.ParticleLayout()
        particle_layout.addParticle(part_coll, 1.0)
        particle_layout.setInterferenceFunction(interference)
        particle_layout.setApproximation(ba.ILayout.SSCA)
        particle_layout.setTotalParticleSurfaceDensity(1)

        # roughness
        r_ptfe = ba.LayerRoughness(2.3*ba.nm, 0.3, 5.0*ba.nm)
        r_hmdso = ba.LayerRoughness(1.1*ba.nm, 0.3, 5.0*ba.nm)

        # layers
        air_layer = ba.Layer(m_air)
        hmdso_layer = ba.Layer(m_HMDSO, self.hmdso_thickness.value)
        hmdso_layer.addLayout(particle_layout)
        ptfe_layer = ba.Layer(m_PTFE, self.ptfe_thickness.value)
        substrate_layer = ba.Layer(m_Si)

        # assembling multilayer
        multi_layer = ba.MultiLayer()
        multi_layer.addLayer(air_layer)
        multi_layer.addLayerWithTopRoughness(hmdso_layer, r_hmdso)
        multi_layer.addLayerWithTopRoughness(ptfe_layer, r_ptfe)
        multi_layer.addLayer(substrate_layer)

        return multi_layer

 

The file FitGALAXIData.py contains parts related to detector initialization, simulation setup, importing the data and finally fit engine setup.

  • The rectangular detector is created to represent PILATUS detector from experiment (line 20).
  • In simulation settings beam is initialized and the detector is assigned to the simulation. Region of interest is assigned at line 40 to simulate only small rectangular window. Additionally, a rectangular mask is added to exclude reflected beam from the analysis (line 42).
  • Real data is loaded from tiff file into histogram representing the detector.
  • run_fitting() function contains the initialization of fitting kernel: creattion of simulation, assignment of fit pair, fit parameters selection (line 62).
  • Two minimizer strategy is selected (Genetic minimizer starts to investigate broad parameter space, Minuit2 minimizer continues to find the local minima.
"""
Real life example: experiment at GALAXY
"""
import matplotlib
from matplotlib import pyplot as plt
import numpy
import bornagain as ba
from SampleBuilder import MySampleBuilder

wavelength = 1.34*ba.angstrom
alpha_i = 0.463*ba.deg

# detector setup as given from instrument responsible
pilatus_npx, pilatus_npy = 981, 1043
pilatus_pixel_size = 0.172  # in mm
detector_distance = 1730.0  # in mm
beam_xpos, beam_ypos = 597.1, 323.4  # in pixels


def create_detector():
    """
    Returns a model of the GALAXY detector
    """
    u0 = beam_xpos*pilatus_pixel_size  # in mm
    v0 = beam_ypos*pilatus_pixel_size  # in mm
    detector = ba.RectangularDetector(pilatus_npx, pilatus_npx*pilatus_pixel_size,
                                      pilatus_npy, pilatus_npy*pilatus_pixel_size)
    detector.setPerpendicularToDirectBeam(detector_distance, u0, v0)
    return detector


def create_simulation():
    """
    Creates and returns GISAS simulation with beam and detector defined
    """
    simulation = ba.GISASSimulation()
    simulation.setDetector(create_detector())
    simulation.setBeamParameters(wavelength, alpha_i, 0.0)
    simulation.setBeamIntensity(1.2e7)
    simulation.setRegionOfInterest(85.0, 70.0, 120.0, 92.)
    # mask on reflected beam
    simulation.addMask(ba.Rectangle(101.9, 82.1, 103.7, 85.2), True)
    # detector resolution function
    # simulation.setDetectorResolutionFunction(
    #   ba.ResolutionFunction2DGaussian(0.5*pilatus_pixel_size,
    #      0.5*pilatus_pixel_size))
    # beam divergence
    # alpha_distr = ba.DistributionGaussian(alpha_i, 0.02*ba.deg)
    # simulation.addParameterDistribution("*/Beam/Alpha", alpha_distr, 5)
    return simulation


def load_real_data(filename="galaxi_data.tif.gz"):
    """
    Fill histogram representing our detector with intensity data from tif file.
    Returns cropped version of it, which represent the area we are interested in.
    """
    hist = ba.IHistogram.createFrom(filename)
    return hist


def run_fitting():
    simulation = create_simulation()
    sample_builder = MySampleBuilder()
    simulation.setSampleBuilder(sample_builder)

    real_data = load_real_data()

    fit_suite = ba.FitSuite()
    draw_observer = ba.DefaultFitObserver(draw_every_nth=10)
    fit_suite.attachObserver(draw_observer)
    fit_suite.initPrint(10)
    fit_suite.addSimulationAndRealData(simulation, real_data)
    print("1.8")

    # setting fitting parameters with starting values
    fit_suite.addFitParameter(
        "*radius", 5.0*ba.nm, ba.AttLimits.limited(4.0, 6.0),
        0.1*ba.nm)
    fit_suite.addFitParameter(
        "*sigma", 0.55, ba.AttLimits.limited(0.2, 0.8), 0.01*ba.nm)
    fit_suite.addFitParameter(
        "*distance", 27.*ba.nm, ba.AttLimits.limited(20, 70),
        0.1*ba.nm)

    use_two_minimizers_strategy = False
    if use_two_minimizers_strategy:
        strategy1 = ba.AdjustMinimizerStrategy("Genetic")
        fit_suite.addFitStrategy(strategy1)

        # Second fit strategy will use another algorithm.
        # It will use best parameters found from previous minimization round.
        strategy2 = ba.AdjustMinimizerStrategy("Minuit2", "Migrad")
        fit_suite.addFitStrategy(strategy2)

    # running fit
    fit_suite.runFit()

    plt.show()

if __name__ == '__main__':
    run_fitting()

Tags: