Skip to content
Snippets Groups Projects
Commit 4af1d504 authored by Pospelov, Gennady's avatar Pospelov, Gennady
Browse files

Beautify fit_galaxi example

parent 70097c85
No related branches found
No related tags found
No related merge requests found
"""
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)
# 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()
"""
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,
ba.RealLimits.limited(0.0, self.hmdso_thickness.value/2.0))
# 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
""" """
Real life example: experiment at GALAXY Fitting experimental data: spherical nanoparticles with size distribution
in 3 layers system (experiment at GALAXI).
""" """
import matplotlib
from matplotlib import pyplot as plt
import numpy
import bornagain as ba import bornagain as ba
from bornagain import nm as nm from bornagain import nm as nm
from SampleBuilder_new import MySampleBuilder from sample_builder import SampleBuilder
wavelength = 1.34*ba.angstrom wavelength = 1.34*ba.angstrom
alpha_i = 0.463*ba.deg alpha_i = 0.463*ba.deg
...@@ -39,17 +37,9 @@ def create_simulation(params): ...@@ -39,17 +37,9 @@ def create_simulation(params):
simulation.setBeamParameters(wavelength, alpha_i, 0.0) simulation.setBeamParameters(wavelength, alpha_i, 0.0)
simulation.setBeamIntensity(1.2e7) simulation.setBeamIntensity(1.2e7)
simulation.setRegionOfInterest(85.0, 70.0, 120.0, 92.) 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) # mask on reflected beam
simulation.addMask(ba.Rectangle(101.9, 82.1, 103.7, 85.2), True)
# detector resolution function sample_builder = SampleBuilder()
# 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)
sample_builder = MySampleBuilder()
sample = sample_builder.create_sample(params) sample = sample_builder.create_sample(params)
simulation.setSample(sample) simulation.setSample(sample)
...@@ -58,11 +48,9 @@ def create_simulation(params): ...@@ -58,11 +48,9 @@ def create_simulation(params):
def load_real_data(filename="galaxi_data.tif.gz"): def load_real_data(filename="galaxi_data.tif.gz"):
""" """
Fill histogram representing our detector with intensity data from tif file. Loads experimental data and returns numpy array.
Returns cropped version of it, which represent the area we are interested in.
""" """
hist = ba.IHistogram.createFrom(filename).array() return ba.IntensityDataIOFactory.readIntensityData(filename).array()
return hist
def run_fitting(): def run_fitting():
......
""" """
3 layers system (substrate, teflon, air). 3 layers system (substrate, teflon, air).
Air layer is populated with spheres with some size distribution. Air layer is populated with spheres with size distribution.
""" """
import bornagain as ba import bornagain as ba
class MySampleBuilder(): class SampleBuilder:
""" """
Builds complex sample from set of parameters.
""" """
def __init__(self): def __init__(self):
# parameters describing the sample """
Init SampleBuilder with default sample parameters.
"""
self.radius = 5.75*ba.nm self.radius = 5.75*ba.nm
self.sigma = 0.4 self.sigma = 0.4
self.distance = 53.6*ba.nm self.distance = 53.6*ba.nm
...@@ -20,13 +22,25 @@ class MySampleBuilder(): ...@@ -20,13 +22,25 @@ class MySampleBuilder():
self.hmdso_thickness = 18.5*ba.nm self.hmdso_thickness = 18.5*ba.nm
def create_sample(self, params): def create_sample(self, params):
"""
Create sample from given set of parameters.
Args:
params: A dictionary containing parameter map.
Returns:
A multi layer.
"""
self.radius = params["radius"] self.radius = params["radius"]
self.sigma = params["sigma"] self.sigma = params["sigma"]
self.distance = params["distance"] self.distance = params["distance"]
return self.multilayer() return self.multilayer()
# constructs the sample for current values of parameters
def multilayer(self): def multilayer(self):
"""
Constructs the sample from current parameter values.
"""
# defining materials # defining materials
m_air = ba.HomogeneousMaterial("Air", 0.0, 0.0) m_air = ba.HomogeneousMaterial("Air", 0.0, 0.0)
m_Si = ba.HomogeneousMaterial("Si", 5.78164736e-6, 1.02294578e-7) m_Si = ba.HomogeneousMaterial("Si", 5.78164736e-6, 1.02294578e-7)
...@@ -38,20 +52,14 @@ class MySampleBuilder(): ...@@ -38,20 +52,14 @@ class MySampleBuilder():
nparticles = 20 nparticles = 20
nfwhm = 2.0 nfwhm = 2.0
sphere_ff = ba.FormFactorFullSphere(self.radius) sphere_ff = ba.FormFactorFullSphere(self.radius)
# sphere_ff = ba.FormFactorTruncatedSphere(
# self.radius, self.radius*1.5)
sphere = ba.Particle(m_Ag, sphere_ff) sphere = ba.Particle(m_Ag, sphere_ff)
position = ba.kvector_t(0*ba.nm, 0*ba.nm, position = ba.kvector_t(0*ba.nm, 0*ba.nm,-1.0*self.hmdso_thickness)
-1.0*self.hmdso_thickness)
sphere.setPosition(position) sphere.setPosition(position)
ln_distr = ba.DistributionLogNormal(self.radius, self.sigma) ln_distr = ba.DistributionLogNormal(self.radius, self.sigma)
par_distr = ba.ParameterDistribution( par_distr = ba.ParameterDistribution(
"/Particle/FullSphere/Radius", ln_distr, nparticles, nfwhm, "/Particle/FullSphere/Radius", ln_distr, nparticles, nfwhm,
ba.RealLimits.limited(0.0, self.hmdso_thickness/2.0)) ba.RealLimits.limited(0.0, self.hmdso_thickness/2.0))
# par_distr = ba.ParameterDistribution(
# "/Particle/TruncatedSphere/Radius", ln_distr, nparticles, nfwhm)
# par_distr.linkParameter("/Particle/TruncatedSphere/Height")
part_coll = ba.ParticleDistribution(sphere, par_distr) part_coll = ba.ParticleDistribution(sphere, par_distr)
# interference function # interference function
...@@ -63,11 +71,11 @@ class MySampleBuilder(): ...@@ -63,11 +71,11 @@ class MySampleBuilder():
interference.setProbabilityDistribution(pdf) interference.setProbabilityDistribution(pdf)
# assembling particle layout # assembling particle layout
particle_layout = ba.ParticleLayout() layout = ba.ParticleLayout()
particle_layout.addParticle(part_coll, 1.0) layout.addParticle(part_coll, 1.0)
particle_layout.setInterferenceFunction(interference) layout.setInterferenceFunction(interference)
particle_layout.setApproximation(ba.ILayout.SSCA) layout.setApproximation(ba.ILayout.SSCA)
particle_layout.setTotalParticleSurfaceDensity(1) layout.setTotalParticleSurfaceDensity(1)
# roughness # roughness
r_ptfe = ba.LayerRoughness(2.3*ba.nm, 0.3, 5.0*ba.nm) r_ptfe = ba.LayerRoughness(2.3*ba.nm, 0.3, 5.0*ba.nm)
...@@ -76,7 +84,7 @@ class MySampleBuilder(): ...@@ -76,7 +84,7 @@ class MySampleBuilder():
# layers # layers
air_layer = ba.Layer(m_air) air_layer = ba.Layer(m_air)
hmdso_layer = ba.Layer(m_HMDSO, self.hmdso_thickness) hmdso_layer = ba.Layer(m_HMDSO, self.hmdso_thickness)
hmdso_layer.addLayout(particle_layout) hmdso_layer.addLayout(layout)
ptfe_layer = ba.Layer(m_PTFE, self.ptfe_thickness) ptfe_layer = ba.Layer(m_PTFE, self.ptfe_thickness)
substrate_layer = ba.Layer(m_Si) substrate_layer = ba.Layer(m_Si)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment