Comparison with PSI reduction software#
This notebook is for validating the results from the essreflectometry packages by comparing them to the ones obtained by the PSI software at the Amor beamline.
[1]:
import orsopy
import matplotlib.pyplot as plt
import scipp as sc
import plopp as pp
from ess import amor
from ess.amor import data # noqa: F401
from ess.reflectometry.types import *
from ess.amor.types import ChopperPhase
Recreate reference files using Jochens Amor data reduction software#
step 1#
%%bash
git clone https://github.com/jochenstahn/amor.git
wget -P raw https://public.esss.dk/groups/scipp/ess/amor/1/amor2023n000{608..614}.hdf
step 2#
%%bash
args='-Y 2023 -d raw -n 614 -m 0.05 -co -7.5 -of ort'
rm raw/614.norm
python amor/neos.py $args -f 614
for i in `seq 608 613`
do
python amor/neos.py $args -f $i -o $i
done
Run essreflectometry on the same files#
[2]:
workflow = amor.AmorWorkflow()
workflow[SampleSize[SampleRun]] = sc.scalar(10.0, unit='mm')
workflow[SampleSize[ReferenceRun]] = sc.scalar(10.0, unit='mm')
workflow[ChopperPhase[SampleRun]] = sc.scalar(-7.5, unit='deg')
# In Jochens workflow:
# * divergence mask: [-0.7, 0.7]
# * note, divergence relative to beam center
workflow[WavelengthBins] = sc.geomspace('wavelength', 2.8, 12.5, 301, unit='angstrom')
workflow[QBins] = sc.geomspace(
dim='Q', start=0.00505, stop=2.93164766e-01, num=391, unit='1/angstrom'
)
workflow[YIndexLimits] = sc.scalar(11, unit=None), sc.scalar(41, unit=None)
workflow[ZIndexLimits] = sc.scalar(80, unit=None), sc.scalar(370, unit=None)
# Chopper phase value in the file is wrong, so we set it manually
workflow[ChopperPhase[ReferenceRun]] = sc.scalar(-7.5, unit='deg')
# The sample rotation value in the file is slightly off, so we set it manually
workflow[SampleRotation[ReferenceRun]] = sc.scalar(0.65, unit='deg')
workflow[Filename[ReferenceRun]] = amor.data.amor_reference_run()
reference_result = workflow.compute(IdealReferenceIntensity)
# Set the result back onto the pipeline to cache it
workflow[IdealReferenceIntensity] = reference_result
/home/runner/work/essreflectometry/essreflectometry/.tox/docs/lib/python3.10/site-packages/scippnexus/nxtransformations.py:222: UserWarning: Failed to convert /entry1/Amor/detector/transformation/height into a transformation: A transformation needs a vector attribute. Falling back to returning underlying value.
transform = Transformation(obj).make_transformation(
/home/runner/work/essreflectometry/essreflectometry/.tox/docs/lib/python3.10/site-packages/scippnexus/nxtransformations.py:222: UserWarning: Failed to convert /entry1/Amor/detector/transformation/rotation into a transformation: A transformation needs a vector attribute. Falling back to returning underlying value.
transform = Transformation(obj).make_transformation(
[3]:
mu = {
'608': 0.8,
'609': 2.2,
'610': 3.6,
'611': 5.0,
'612': 0.6,
'613': 0.6,
}
results = sc.DataGroup({'ess': sc.DataGroup(), 'psi': sc.DataGroup()})
# ESS results
for key, angle in mu.items():
print(key, '... ', end='')
workflow[SampleRotation[SampleRun]] = sc.scalar(angle + 0.05, unit='deg')
workflow[Filename[SampleRun]] = amor.data.amor_sample_run(key)
da = workflow.compute(NormalizedIofQ).hist()
da.coords['Q'] = sc.midpoints(da.coords['Q'], dim='Q')
results['ess'][key] = da
print('Done!')
608 ... Done!
609 ... Done!
610 ... Done!
611 ...
Downloading file 'amor2023n000612.hdf' from 'https://public.esss.dk/groups/scipp/ess/amor/2/amor2023n000612.hdf' to '/home/runner/.cache/ess/amor/2'.
Done!
612 ...
Downloading file 'amor2023n000613.hdf' from 'https://public.esss.dk/groups/scipp/ess/amor/2/amor2023n000613.hdf' to '/home/runner/.cache/ess/amor/2'.
Done!
613 ... Done!
[4]:
# PSI results
for key in mu:
d = orsopy.fileio.load_orso(amor.data.amor_psi_software_result(key))[0].data
y = sc.array(dims=['Q'], values=d[:, 1], variances=d[:, 2] ** 2)
results['psi'][key] = sc.DataArray(
data=y, coords={'Q': sc.array(dims=['Q'], values=d[:, 0], unit='1/angstrom')}
)
Downloading file '608.Rqz.ort' from 'https://public.esss.dk/groups/scipp/ess/amor/2/608.Rqz.ort' to '/home/runner/.cache/ess/amor/2'.
Downloading file '609.Rqz.ort' from 'https://public.esss.dk/groups/scipp/ess/amor/2/609.Rqz.ort' to '/home/runner/.cache/ess/amor/2'.
Downloading file '610.Rqz.ort' from 'https://public.esss.dk/groups/scipp/ess/amor/2/610.Rqz.ort' to '/home/runner/.cache/ess/amor/2'.
Downloading file '611.Rqz.ort' from 'https://public.esss.dk/groups/scipp/ess/amor/2/611.Rqz.ort' to '/home/runner/.cache/ess/amor/2'.
Downloading file '612.Rqz.ort' from 'https://public.esss.dk/groups/scipp/ess/amor/2/612.Rqz.ort' to '/home/runner/.cache/ess/amor/2'.
Downloading file '613.Rqz.ort' from 'https://public.esss.dk/groups/scipp/ess/amor/2/613.Rqz.ort' to '/home/runner/.cache/ess/amor/2'.
Plot reflectivity curve comparison#
[5]:
fig, axs = plt.subplots(3, 2, figsize=(12, 12))
for ax, (key, angle) in zip(axs.ravel(), mu.items(), strict=True):
pp.plot(
{
facility: results[facility][key][results[facility][key].data > 0.0]
for facility in results
},
ax=ax,
norm='log',
ls='solid',
marker=None,
vmin=1.0e-4,
title=f'Angle: {angle} deg',
)
fig.tight_layout()
/home/runner/work/essreflectometry/essreflectometry/.tox/docs/lib/python3.10/site-packages/matplotlib/axes/_axes.py:3749: RuntimeWarning: invalid value encountered in add
low, high = dep + np.vstack([-(1 - lolims), 1 - uplims]) * err
/home/runner/work/essreflectometry/essreflectometry/.tox/docs/lib/python3.10/site-packages/plopp/backends/matplotlib/line.py:173: RuntimeWarning: invalid value encountered in subtract
arr2 = np.array([y - e, y + e]).T.flatten()
Overlay curves from same sample at different angle, should be on top of each other#
[6]:
fig, axs = plt.subplots(1, 2, figsize=(12, 4))
for ax, (facility, res) in zip(axs.ravel(), results.items(), strict=True):
pp.plot(
dict(list(res.items())[:4]), # First 4 curves
ax=ax,
norm='log',
ls='solid',
marker=None,
vmin=1.0e-4,
title=facility,
)
