Find peaks

This example simulates GISAS from a one-dimensional grating of long boxes and then applies FindPeaks to the resulting intensity map. The detected peak coordinates are overlaid on the detector image, demonstrating how peak finding can be added as a post-processing step after a scattering simulation.

The simulation also installs a terminal progress monitor for interactive runs. It skips the monitor when __no_terminal__ is present, so the same script can run in batch environments without terminal output.

Result

Find peaks result

Sample

Find peaks sample

Python script

 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
#!/usr/bin/env python3
# /// script
# requires-python = ">=3.10"
# dependencies = ["bornagain>=25,<26"]
# ///
"""
Simulation of grating using very long boxes and 1D lattice.
Monte-carlo integration is used to get rid of
large-particle form factor oscillations.
"""
import bornagain as ba
ba.require_versions("bornagain>=25,<26")
from bornagain import angstrom, deg, micrometer, nm


def get_sample(lattice_rotation_angle=0*deg):
    """
    A sample with a grating on a substrate.
    lattice_rotation_angle = 0 - beam parallel to grating lines
    lattice_rotation_angle = 90*deg - beam perpendicular to grating lines
    """
    # Material
    si_color = (0.30, 0.62, 0.86)
    si_mat = ba.RefractiveMaterial("Si", si_color, 5.7816e-6, 1.0229e-7)

    box_length, box_width, box_height = 50*micrometer, 70*nm, 50*nm
    lattice_length = 150*nm

    box_ff = ba.LongBoxLorentz(box_length, box_width, box_height)
    box = ba.Particle(si_mat, box_ff)
    box.rotate(ba.RotationZ(lattice_rotation_angle))

    # collection of particles
    layout = ba.Crystal1D(box, lattice_length,
                          90*deg - lattice_rotation_angle,
                          1/box_length)
    profile = ba.Profile1DGauss(450)
    layout.setDecayFunction(profile)

    sigma, hurst, corrLength = 5*nm, 0.5, 10*nm
    autocorr = ba.SelfAffineFractalModel(sigma, hurst, corrLength)
    transient = ba.TanhTransient()
    roughness = ba.Roughness(autocorr, transient)

    # assembling the sample
    vacuum_layer = ba.Layer(ba.Vacuum())
    vacuum_layer.deposit2D(layout)
    substrate_layer = ba.Layer(si_mat, roughness)

    sample = ba.Sample()
    sample.addLayer(vacuum_layer)
    sample.addLayer(substrate_layer)
    return sample


def get_simulation(sample):
    beam = ba.Beam(1e8, 1.34*angstrom, 0.4*deg)
    n = 401
    detector = ba.SphericalDetector(n, -0.5*deg, 0.5*deg, n, 0, 0.5*deg)
    simulation = ba.ScatteringSimulation(beam, sample, detector)
    simulation.options().setMonteCarloIntegration(True, 100, seed=0)
    return simulation


if __name__ == '__main__':
    sample = get_sample()
    simulation = get_simulation(sample)
    if not "__no_terminal__" in globals():
        simulation.setTerminalProgressMonitor()
    result = simulation.simulate()
    peaks = ba.FindPeaks(result, 2, "nomarkov", 0.001)
    xpeaks = [peak[0] for peak in peaks]
    ypeaks = [peak[1] for peak in peaks]
    print(peaks)

    ba.plt.plot(xpeaks,
             ypeaks,
             marker='x',
             linestyle='none',
             color='white',
             markersize=10)

    ba.showSample3D(sample, sample_size=1500*nm, seed=0)
    ba.plot_datafield(result, unit_aspect=1)
    ba.plt.show()
auto/Examples/gisas/methods/FindPeaks.py