Accessing simulation results

This is an extended example for the further treatment of simulation results: accessing the results, plotting, cropping, slicing and exporting. This serves as a supporting example to the Accessing simulation results tutorial.

  • The standard Cylinders in DWBA sample is used for running the simulation.
  • The simulation results are retrieved as a Histogram2D object and then processed in various functions to achieve a resulting image.

Intensity images

  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
#!/usr/bin/env python3
"""
Extended example for simulation results treatment (cropping, slicing, exporting)
"""
import math
import random
import bornagain as ba
from bornagain import angstrom, ba_plot as bp, deg, nm, std_samples
from matplotlib import pyplot as plt
import datetime

def get_sample():
    return std_samples.cylinders()


def get_simulation(sample):
    """
    A GISAXS simulation with beam and detector defined.
    """
    beam = ba.Beam(1e5, 1*angstrom, 0.2*deg)
    n = 200
    det = ba.SphericalDetector(n, -2*deg, 2*deg, n, 0, 2*deg)
    simulation = ba.ScatteringSimulation(beam, sample, det)
    return simulation


def get_noisy_image(field):
    """
    Returns clone of input field filled with additional noise
    """
    result = field.clone()
    noise_factor = 2.0
    for i in range(0, result.size()):
        amplitude = field.valAt(i)
        sigma = noise_factor*math.sqrt(amplitude)
        noisy_amplitude = random.gauss(amplitude, sigma)
        result.setAt(i, noisy_amplitude)
    return result


def plot_histogram(field, **kwargs):
    bp.plot_histogram(field,
                      xlabel=r'$\varphi_f ^{\circ}$',
                      ylabel=r'$\alpha_{\rm f} ^{\circ}$',
                      zlabel="",
                      **kwargs)


def plot_slices(noisy):
    """
    Plot 1D slices along y-axis at certain x-axis values.
    """
    plt.yscale('log')

    # projection along Y, slice at fixed x-value
    proj1 = noisy.yProjection(0)
    plt.plot(proj1.axis(0).binCenters(),
             proj1.flatVector(),
             label=r'$\varphi=0.0^{\circ}$')

    # projection along Y, slice at fixed x-value
    proj2 = noisy.yProjection(0.5)  # slice at fixed value
    plt.plot(proj2.axis(0).binCenters(),
             proj2.flatVector(),
             label=r'$\varphi=0.5^{\circ}$')

    # projection along Y for all X values between [xlow, xup], averaged
    proj3 = noisy.yProjection(0.41, 0.59)
    plt.plot(proj3.axis(0).binCenters(),
             proj3.flatVector(),
             label=r'$<\varphi>=0.5^{\circ}$')

    plt.xlim(proj1.axis(0).min(), proj1.axis(0).max())
    plt.ylim(proj1.maxVal()*3e-6, proj1.maxVal()*3)
    plt.xlabel(r'$\alpha_{\rm f} ^{\circ}$', fontsize=16)
    plt.legend(loc='upper right')


def plot(field):
    """
    Demonstrates modified data plots.
    """
    plt.subplots(2, 2, figsize=(12.80, 10.24))

    plt.subplot(2, 2, 1)
    bp.plot_histogram(field)
    plt.title("Intensity as heatmap")

    plt.subplot(2, 2, 2)
    crop = field.crop(-1, 0.5, 1, 1)
    bp.plot_histogram(crop)
    plt.title("Cropping")

    plt.subplot(2, 2, 3)
    noisy = get_noisy_image(field)
    reldiff = ba.relativeDifferenceField(noisy, field).npArray()
    bp.plot_array(reldiff, intensity_min=1e-03, intensity_max=10)
    plt.title("Relative difference")

    plt.subplot(2, 2, 4)
    plot_slices(noisy)
    plt.title("Various slicing of 2D into 1D")

    plt.tight_layout()

if __name__ == '__main__':
    sample = get_sample()
    simulation = get_simulation(sample)
    result = simulation.simulate()

    if bp.datfile:
        ba.writeDatafield(result, bp.datfile + ".int.gz")
        # Other supported extensions are .tif and .txt.
        # Besides compression .gz, we support .bz2, and uncompressed.

    field = result.extracted_field()
    plot(field)
    bp.show_or_export()
auto/Examples/varia/AccessingSimulationResults.py