Reduction of spin flip data from McStas simulation#

This notebook contains an example of polarized reflectometry data reduction. The dataset comes from a McStas simulation of the ESTIA instrument.

Four samples were simulated:

supermirror - a perfect nonmagnetic supermirror
magnetic_supermirror - a magnetic supermirror with high reflectivity for one spin state and low reflectivity for the other
magnetic_supermirror_2 - a magnetic supermirror with slightly different reflectivity curve than magnetic_supermirror
spin_flip_sample - a magnetic supermirror that also causes some incident neutrons to spin flip with 10% probability

Each sample was measured for four different flipper settings offoff, offon, onoff, onon, corresponding to the polarizer and analyzer flipper setting.

In the data reduction procedure supermirror and magnetic_supermirror_2 will be used to calibrate the parameters of the instrument.

We start by importing packages and creating the workflow:

[1]:

#%matplotlib ipympl
import scipp as sc

from ess.reflectometry.types import *
from ess.reflectometry.supermirror import CriticalEdge

from ess.estia.types import *
from ess.estia import EstiaMcStasWorkflow
from ess.estia.data import (
    estia_mcstas_spin_flip_example,
    estia_mcstas_spin_flip_example_groundtruth,
    estia_mcstas_spin_flip_example_download_all_to_cache,
    estia_tof_lookup_table,
)
from ess.reflectometry.figures import wavelength_z_figure
from ess.estia.mcstas import mcstas_wavelength_coordinate_transformation_graph

from ess.polarization import (  # noqa: F401
    Up, Down, ReducedSampleDataBySpinChannel, SupermirrorEfficiencyFunction, Polarizer, Analyzer, PolarizationAnalysisWorkflow,
    SupermirrorWorkflow, SecondDegreePolynomialEfficiency, EfficiencyLookupTable
)
from ess.polarization.types import TotalPolarizationCorrectedData

wf = EstiaMcStasWorkflow()
wf[YIndexLimits]  = sc.scalar(35), sc.scalar(64)
wf[ZIndexLimits] = sc.scalar(600), sc.scalar(930)
wf[BeamDivergenceLimits] = sc.scalar(-0.75, unit='deg'), sc.scalar(0.75, unit='deg')
wf[WavelengthBins] = sc.geomspace('wavelength', 4., 12, 700, unit='angstrom')
wf[QBins] = sc.linspace('Q', 0.001, 0.1, 200, unit='1/angstrom')


wf.insert( mcstas_wavelength_coordinate_transformation_graph )
wf[TimeOfFlightLookupTableFilename] = estia_tof_lookup_table()

# Reference sample is perfect supermirror with reflectivity = 1 everywhere
wf[CriticalEdge] = sc.scalar(float('inf'), unit='1/angstrom')

# There is no proton current data in the McStas files, here we just add some fake proton current
# data to make the workflow run.
wf[ProtonCurrent[SampleRun]] = sc.DataArray(
    sc.array(dims=('time',), values=[]),
    coords={'time': sc.array(dims=('time',), values=[], unit='s')})
wf[ProtonCurrent[ReferenceRun]] = sc.DataArray(
    sc.array(dims=('time',), values=[]),
    coords={'time': sc.array(dims=('time',), values=[], unit='s')})

Download the data: (might take ~2 minutes depending on your internet connection)

[2]:

%%time
estia_mcstas_spin_flip_example_download_all_to_cache()

Downloading file 'spin_flip_example/supermirror_offoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/supermirror_offoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/supermirror_onoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/supermirror_onoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/supermirror_onon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/supermirror_onon.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_offoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_offoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_offon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_offon.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_onon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_onon.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/supermirror_offon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/supermirror_offon.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_2_offoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_2_offoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_onoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_onoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_2_onon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_2_onon.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_2_offon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_2_offon.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/spin_flip_sample_onoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/spin_flip_sample_onoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/spin_flip_sample_offoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/spin_flip_sample_offoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/spin_flip_sample_offon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/spin_flip_sample_offon.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/magnetic_supermirror_2_onoff.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/magnetic_supermirror_2_onoff.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/spin_flip_sample_onon.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/spin_flip_sample_onon.h5' to '/home/runner/.cache/ess/estia'.

CPU times: user 3min 39s, sys: 3min 16s, total: 6min 55s
Wall time: 4min 34s

Reducing the data#

First each dataset is loaded and reduced separately. The datasets are reduced “as references” or “as samples” depending on how they are supposed to be used. supermirror and magnetic_supermirror_2 are used as references to calibrate the instrument, and spin_flip_sample and magnetic_supermirror are reduced as samples.

[3]:

references = {}
for sample in (
    'supermirror',
    'magnetic_supermirror_2',
):
    references[sample] = []
    for flipper_setting in ('offoff', 'offon', 'onoff', 'onon'):
        w = wf.copy()
        w[RawDetector[ReferenceRun]] = sc.io.load_hdf5(estia_mcstas_spin_flip_example(sample, flipper_setting))
        references[sample].append(w.compute(ReducedReference))

        # We need to unalign all coords of the references to use
        # them in the calibration procedure
        for c in references[sample][-1].coords:
            references[sample][-1].coords.set_aligned(c, False)


samples = {}
for sample in (
    'spin_flip_sample',
    'magnetic_supermirror'
):
    samples[sample] = []
    for flipper_setting in ('offoff', 'offon', 'onoff', 'onon'):
        w = wf.copy()
        w[RawDetector[SampleRun]] = sc.io.load_hdf5(estia_mcstas_spin_flip_example(sample, flipper_setting))
        samples[sample].append(w.compute(ReducibleData[SampleRun]))

        # We need to unalign all coords
        for c in samples[sample][-1].coords:
            samples[sample][-1].coords.set_aligned(c, False)

Here we load the ground truth reflectivity curves for the up respectively down spin state. Those will be used to compare to the computed reflectivity curves.

[4]:

Rdown = sc.io.load_hdf5(estia_mcstas_spin_flip_example_groundtruth('down'))
Rup = sc.io.load_hdf5(estia_mcstas_spin_flip_example_groundtruth('up'))

Downloading file 'spin_flip_example/ground_truth_spin_down_reflectivity.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/ground_truth_spin_down_reflectivity.h5' to '/home/runner/.cache/ess/estia'.
Downloading file 'spin_flip_example/ground_truth_spin_up_reflectivity.h5' from 'https://public.esss.dk/groups/scipp/ess/estia/1/spin_flip_example/ground_truth_spin_up_reflectivity.h5' to '/home/runner/.cache/ess/estia'.

To calibrate the workflow we estimate the polarization efficiency of the polarizer and the analyzer respectively.

The calibration is performed on the wavelength-dependent intensity obtained by summing over a slice of the detector. We don’t sum over the entire detector because the signal is best in the center region. In the lower part of the detector the reflectivity of the magnetic reference increases quickly to be the same as the reflectivity of the non magnetic supermirror reference. In that region we cannot distinguish between the spin states, so doing the calibration there is useless. In the top part of the detector the signal is lower and there is more noise. We don’t want to make the region of the detecor we consider too small because then we will have too little signal.

The trade of is that we only consider the region of the detector highlighted in the below figure:

[5]:

wavelength_z_figure(references['magnetic_supermirror_2'][0].assign_masks(
    bad_region=sc.where(
        (references['magnetic_supermirror_2'][0].coords['z_index'] < sc.scalar(21 * 32)) |
        (references['magnetic_supermirror_2'][0].coords['z_index'] >= sc.scalar(25 * 32)),
        sc.scalar(True),
        sc.scalar(False)
    )
))

[5]:

../../_images/user-guide_estia_simulated-spin-flip-sample_10_0.svg

The intensity is accumulated over each blade of the detector, corresponding to a relatively thin interval of scattering angle. In such a thin interval of scattering angle the reflectivity of the magnetic supermirror is roughly a function of wavelength, and this is what we need for the calibration to work well.

From the procedure described above we obtain separate calibration results for each blade of the detector. The separate calibration values are combined using a weighted average where each calibration value is weighted by the inverse of its variance.

[6]:

from ess.estia.calibration import PolarizationCalibrationParameters


calibration = PolarizationCalibrationParameters.from_reference_measurements(
    [r['blade', 21:25].sum('wire').rebin(wavelength=sc.linspace('wavelength', 4, 10.7, 50, unit='angstrom')) for r in references['supermirror']],
    [r['blade', 21:25].sum('wire').rebin(wavelength=sc.linspace('wavelength', 4, 10.7, 50, unit='angstrom')) for r in references['magnetic_supermirror_2']],
)


def weighted_mean(p, dim):
    return ((p / sc.variances(p)).nanmean(dim) / ((1 / sc.variances(p)).nanmean(dim)))


# Assuming the flipper efficiency is wavelength independent (as is the case in this simulation) we can ignore the other parameters
polarizer_efficiency = weighted_mean(calibration.Pp, 'blade')
analyzer_efficiency = weighted_mean(calibration.Ap, 'blade')

[7]:

polarizer_efficiency.plot(title='Polarizer efficiency') + analyzer_efficiency.plot(title='Analyzer efficiency')

[7]:

../../_images/user-guide_estia_simulated-spin-flip-sample_13_0.svg

Finally the calibration curves are used to fit quadratic polynomials that approximate the polarization efficiency as functions of wavelength.

[8]:

def efficiency_model(wavelength, a, b, c):
    '''The polarization efficiency of the analyzer and the polarizer
    is modeled as a quadratic polynomial'''
    return a * wavelength**2 + b * wavelength + c


def fit_model_to_efficiency(efficiency):
    da = efficiency.copy()
    da.coords['wavelength'] = sc.midpoints(da.coords['wavelength'])
    par, _ = sc.curve_fit(
        ['wavelength'],
        efficiency_model,
        da,
        p0={
            'a': sc.scalar(1., unit='1/angstrom**2'),
            'b': sc.scalar(1., unit='1/angstrom'),
            'c': sc.scalar(1., unit='dimensionless')
        }
    )
    return {k: sc.values(p.data) for k, p in par.items()}


sc.plot({
    'fit': sc.DataArray(
        efficiency_model(
            polarizer_efficiency.coords['wavelength'],
            **fit_model_to_efficiency(polarizer_efficiency)
        ),
        coords=polarizer_efficiency.coords
    ),
    'efficiency': polarizer_efficiency
}, title='polarizer') + sc.plot({
    'fit': sc.DataArray(
        efficiency_model(
            analyzer_efficiency.coords['wavelength'],
            **fit_model_to_efficiency(analyzer_efficiency)
        ),
        coords=analyzer_efficiency.coords
    ),
    'efficiency': analyzer_efficiency
}, title='analyzer')

[8]:

../../_images/user-guide_estia_simulated-spin-flip-sample_15_0.svg

In the polarization correction workflow we can use either the obtained calibration values directly (as a lookup table in wavelength) or we can use the polynomial curves to compute polarization efficiencies as functions of wavelength.

[9]:

pwf = PolarizationAnalysisWorkflow(polarizer_workflow=SupermirrorWorkflow(), analyzer_workflow=SupermirrorWorkflow())

# To do the calibration using a lookup table - uncomment these lines
#pwf[SupermirrorEfficiencyFunction[Polarizer]] = EfficiencyLookupTable(polarizer_efficiency)
#pwf[SupermirrorEfficiencyFunction[Analyzer]] = EfficiencyLookupTable(analyzer_efficiency)

# To do the calibration using the polynomial approximation - uncomment these lines
pwf[SupermirrorEfficiencyFunction[Polarizer]] = SecondDegreePolynomialEfficiency(**fit_model_to_efficiency(polarizer_efficiency))
pwf[SupermirrorEfficiencyFunction[Analyzer]] = SecondDegreePolynomialEfficiency(**fit_model_to_efficiency(analyzer_efficiency))


for i, pols in enumerate((Down, Up)):
    for j, anas in enumerate((Down, Up)):
        sam = references['supermirror'][2*i + j]
        pwf[ReducedSampleDataBySpinChannel[pols, anas]] = sam.assign_coords(wavelength=sc.midpoints(sam.coords['wavelength']))


res = pwf.compute(TotalPolarizationCorrectedData)
I0 = (res.upup + res.downdown) / 2
mask = sc.isnan(I0).sum(('blade', 'wire')) > sc.scalar(0, unit=None)
wf[ReducedReference] = I0.assign_masks(nopolcal=mask.data).assign_coords(wavelength=wf.compute(WavelengthBins))

# Required to read sample rotation / similar parameters associated with the reference measurement
wf[RawDetector] = sc.io.load_hdf5(estia_mcstas_spin_flip_example('supermirror', 'offoff'))

sample_name = 'spin_flip_sample'

for i, pols in enumerate((Down, Up)):
    for j, anas in enumerate((Down, Up)):
        w = wf.copy()
        w[ReducibleData[SampleRun]] = samples[sample_name][2*i + j]
        sam = w.compute(ReflectivityOverQ).bin(wavelength=wf.compute(WavelengthBins)).assign_masks(nopolcal=mask.data)
        pwf[ReducedSampleDataBySpinChannel[pols, anas]] = sam.assign_coords(wavelength=sc.midpoints(sam.coords['wavelength']))


spin_flip_sample_reflectivity = pwf.compute(TotalPolarizationCorrectedData)


sample_name = 'magnetic_supermirror'

for i, pols in enumerate((Down, Up)):
    for j, anas in enumerate((Down, Up)):
        w = wf.copy()
        w[ReducibleData[SampleRun]] = samples[sample_name][2*i + j]
        sam = w.compute(ReflectivityOverQ).bin(wavelength=wf.compute(WavelengthBins)).assign_masks(nopolcal=mask.data)
        pwf[ReducedSampleDataBySpinChannel[pols, anas]] = sam.assign_coords(wavelength=sc.midpoints(sam.coords['wavelength']))


magnetic_supermirror_reflectivity = pwf.compute(TotalPolarizationCorrectedData)


(sc.plot(
    {'Rdd': spin_flip_sample_reflectivity.downdown.sum(('wavelength',)),
     'Ruu': spin_flip_sample_reflectivity.upup.sum(('wavelength',)),
     'Rdu': spin_flip_sample_reflectivity.downup.sum(('wavelength',)),
     'Rud': spin_flip_sample_reflectivity.updown.sum(('wavelength',)),
     'True Rdd': Rdown * 0.9,
     'True Ruu': Rup * 0.9,
     'True Rdu': Rdown * 0.1,
    }, norm='log', title='Reflectivity of spin flip sample', vmin=1e-4) +
sc.plot(
    {'Rdd': magnetic_supermirror_reflectivity.downdown.sum(('wavelength',)),
     'Ruu': magnetic_supermirror_reflectivity.upup.sum(('wavelength',)),
     'Rdu': magnetic_supermirror_reflectivity.downup.sum(('wavelength',)),
     'Rud': magnetic_supermirror_reflectivity.updown.sum(('wavelength',)),
     'True Rdd': Rdown,
     'True Ru': Rup,
     'True Rdu': Rdown * 0.0,
    }, norm='log', title='Reflectivity of magnetic supermirror sample', vmin=1e-4))

[9]:

../../_images/user-guide_estia_simulated-spin-flip-sample_17_0.svg

This Page

Reduction of spin flip data from McStas simulation#

Reducing the data#