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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297

#!/usr/bin/env python3
"""
Real life example:
Fitting data from an X-ray reflectometer

The sample is composed of a thin
silver nano-particle layer on a silicon
substrate. The substrate is covered with
SiO2 layer. The nano-particle layer has negligible density
and does not considerably affect
the observed reflectometry picture.

The following parameters of the experiment
are fitted:
1. beam intensity
2. footprint correction factor
3. beam angular divergence
4. Material concentration in the SiO2 layer
5. Thickness of SiO2 layer
6. Sample roughness

Fitting is done in two steps:
First the whole range of experimental data is fitted,
then the data related to the instrument is fixed and
on the second step only the right-hand part of
experimental data (i.e. the part of the reflectometry curve
associated with bigger incident angles)
is concerned for fitting. At the second
stage only the sample parameters are fitted,
since only they affect the shape of the reflectometry
curve at bigger incident angles.
"""

import os
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import differential_evolution
import bornagain as ba
from bornagain import ba_plot


def get_real_data(filename):
    """
    Loads real data files and merges them once.
    Returns a Nx3 array (N - the number of experimental data entries)
    with first column being coordinates,
    second one being values,
    and the third one being weights to restore intensity values from experiment
    """

    real_data = np.loadtxt(filename, usecols=(0, 1, 3), skiprows=1)

    # translating axis values from double incident angle (degs)
    # to incident angle (radians)
    real_data[:, 0] *= np.pi/360

    global expdata
    expdata = real_data.copy()


def get_real_data_axis(start, end):
    """
    Get axis coordinates of the experimental data
    :param start: first bin to extract
    :param end: last bin to extract
    :return: 1D array with axis coordinates
    """
    return expdata[start:end, 0]


def get_real_data_values(start, end):
    """
    Get experimental data values as a 1D array
    :param start: first bin to extract
    :param end: last bin to extract
    :return: 1D array with experimental data values
    """
    return expdata[start:end, 1]


def get_weights(start, end):
    """
    Get weights to restore genuine intensity of experimental instrument
    :param start: first bin to extract
    :param end: last bin to extract
    :return: 1D array with weights to restore beam intensity
    """
    return expdata[start:end, 2]


def create_simulation(sample, arg_dict, bin_start, bin_end):
    """
    Creates and returns specular simulation
    """
    wavelength = 1.54*ba.angstrom
    alpha_distr = ba.DistributionGaussian(0, arg_dict["divergence"])
    footprint = ba.FootprintGauss(arg_dict["footprint_factor"])

    scan = ba.AlphaScan(wavelength, get_real_data_axis(bin_start, bin_end))
    scan.setAngleDistribution(alpha_distr)
    scan.setFootprintFactor(footprint)

    simulation = ba.SpecularSimulation(scan, sample)
    return simulation


def buildSample(arg_dict):
    """
    Creates sample and returns it
    """
    # defining materials
    m_vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
    m_si_o2 = ba.RefractiveMaterial(
        "SiO2", 8.57040868e-06*arg_dict["concentration"],
        1.11016654e-07*arg_dict["concentration"])
    m_si = ba.RefractiveMaterial("Si", 7.57211137e-06, 1.72728178e-07)

    # roughness
    r_si = ba.LayerRoughness(arg_dict["roughness"], 0, 0)

    # layers
    vacuum_layer = ba.Layer(m_vacuum)
    oxide_layer = ba.Layer(m_si_o2, arg_dict["thickness"])
    substrate_layer = ba.Layer(m_si)

    # assembling sample
    sample = ba.MultiLayer()
    sample.addLayer(vacuum_layer)
    sample.addLayerWithTopRoughness(oxide_layer, r_si)
    sample.addLayerWithTopRoughness(substrate_layer, r_si)

    return sample


def run_simulation(arg_dict, bin_start=0, bin_end=-1):
    """
    Runs simulation and returns its result
    """

    sample = buildSample(arg_dict)
    simulation = create_simulation(sample, arg_dict, bin_start, bin_end)

    result = simulation.simulate()
    result.data_field().scale(arg_dict["intensity"])

    return result


def chi_2(real_data, sim_data, weights):
    """
    Computes chi_2 metrics and returns its value
    """
    sim_data_upsc = np.multiply(weights, sim_data)
    sim_data_upsc[sim_data_upsc == 0] = 1e-30
    real_data_upsc = np.multiply(weights, real_data)
    diff = real_data_upsc - sim_data_upsc
    return np.sum(np.divide(np.multiply(diff, diff), sim_data_upsc))


def create_par_dict(*arg):
    """
    Creates a dictionary with parameter names and values
    and returns it
    """
    return {
        'intensity': arg[0],
        'footprint_factor': arg[1],
        'divergence': arg[2],
        'concentration': arg[3],
        'thickness': arg[4],
        'roughness': arg[5]
    }


def objective_primary(args):
    """
    Objective function for preliminary stage of optimization
    """

    bin_start = 15  # first bin in the experimental data to calculate
    bin_end = -1  # last bin in the experimental data to calculate
    arg_dict = create_par_dict(*args)

    sim_result = run_simulation(arg_dict, bin_start, bin_end)
    sim_data = sim_result.array()
    return chi_2(get_real_data_values(bin_start, bin_end), sim_data,
                 get_weights(bin_start, bin_end))


def objective_fine(args, intensity, footprint_factor, divergence):
    """
    Objective function for tuning the right-hand side of experimental data
    """

    bin_start = 404  # first bin in the experimental data to calculate
    bin_end = -1  # last bin in the experimental data to calculate
    arg_dict = create_par_dict(intensity, footprint_factor, divergence,
                               *args)

    sim_result = run_simulation(arg_dict, bin_start, bin_end)
    sim_data = sim_result.array()
    return chi_2(get_real_data_values(bin_start, bin_end), sim_data,
                 get_weights(bin_start, bin_end))


def run_fitting():
    """
    Runs fitting and returns its result
    """
    # running preliminary optimization on the total range of experimental data.
    bounds = [
        (1e6, 1e8),  # beam intensity
        (0, 0.1),  # beam-to-sample width ratio
        (0, 0.08*ba.deg),  # beam_divergence
        (0, 1),  # oxide_concentration
        (0, 2*ba.nm),  # oxide_thickness
        (0, 2*ba.nm)
    ]  # roughness

    print("Start preliminary fitting of experimental data:\n")

    preliminary_result = differential_evolution(objective_primary,
                                                bounds,
                                                maxiter=20,
                                                popsize=60,
                                                mutation=(0.5, 1.5),
                                                disp=True,
                                                tol=1e-5)

    bounds = [
        (0, 1),  # oxide_concentration
        (0, 2*ba.nm),  # oxide_thickness
        (0, 2*ba.nm)
    ]  # roughness

    fixed_args = (
        preliminary_result.x[0],  # beam intensity
        preliminary_result.x[1],  # beam-to-sample width ratio
        preliminary_result.x[2]  # beam divergence
    )

    print(
        "\nStart fitting big incident angle part of experimental data:\n")

    fine_tuning_result = differential_evolution(objective_fine,
                                                bounds,
                                                fixed_args,
                                                maxiter=20,
                                                popsize=40,
                                                mutation=(0.5, 1.5),
                                                disp=True,
                                                tol=1e-5)

    result = create_par_dict(*fixed_args, *fine_tuning_result.x)
    print("\nFitting result:")
    print(result, "\n")

    return result


def plot_result(sim_result, ref_result, bin_start=0, bin_end=-1):
    """
    Plots the graphs of obtained simulation data
    """
    sim_data = sim_result.array()
    ref_data = ref_result.array()

    plt.yscale('log')
    plt.plot(
        get_real_data_axis(bin_start, bin_end)*180/np.pi,
        get_real_data_values(bin_start, bin_end), sim_result.axis(),
        sim_data, ref_result.axis(), ref_data)

    xlabel = ba_plot.get_axes_labels(sim_result, ba.Coords_UNDEFINED)[0]
    ylabel = "Intensity"
    plt.xlabel(xlabel, fontsize=16)
    plt.ylabel(ylabel, fontsize=16)
    plt.legend(['Experimental data', 'Simulation', 'Reference'],
               loc='upper right',
               fontsize=16)

    plt.show()


if __name__ == '__main__':
    datadir = os.getenv('BA_EXAMPLE_DATA_DIR', '')
    data_fname = os.path.join(datadir, "mg6a_Merged.txt.gz")
    get_real_data(data_fname)
    fit_data = run_fitting()
    ref_data = create_par_dict(
        3.78271438e6,  # beam intensity
        9.58009763e-04,  # beam-to-sample width ratio
        2.30471294e-04,  # beam angular divergence
        0.58721753,  # oxide concentration
        1.25559347,  # oxide thickness
        0.19281863)  # roughness
    plot_result(run_simulation(fit_data), run_simulation(ref_data))