Zoom Polarization Analysis

Introduction

%matplotlib widget
import sciline
import scipp as sc
from ess import polarization as pol
from ess import sans
from ess import isissans as isis
from ess.sans.types import *

sans_workflow = isis.zoom.ZoomWorkflow()
sans_workflow = sans.with_pixel_mask_filenames(sans_workflow, [])  # no masks

from pathlib import Path

data_folder = Path('zoom_polarized_data')
# Runs with analyzer at 4 different times
cell_runs = [str(data_folder / f'ZOOM00022{run}.nxs') for run in [710, 712, 714, 716]]
empty_run = data_folder / 'ZOOM00034787.nxs'
depolarized_run = data_folder / 'ZOOM00022718.nxs'
cell_runs = cell_runs + [depolarized_run]

Setup SANS workflow for computing transmission fraction

from ess.polarization.zoom import ZoomTransmissionFractionWorkflow

sans_workflow = ZoomTransmissionFractionWorkflow(cell_runs)
sans_workflow[Filename[EmptyBeamRun]] = str(empty_run)
sans_workflow[WavelengthBins] = sc.geomspace(
    'wavelength', start=1.75, stop=16.5, num=141, unit='Å'
)
sans_workflow[isis.mantidio.Period] = 0
sans_workflow[isis.mantidio.CalibrationWorkspace] = None

Inspect data for one of the runs with analyzer cell

loaded = sciline.get_mapped_node_names(
    sans_workflow, isis.mantidio.LoadedFileContents[TransmissionRun[SampleRun]]
)
# Load only first run
first_run = sans_workflow.compute(loaded[0])
sc.DataGroup(sc.collapse(first_run['data'], keep='tof')).plot()

first_run

We can load the combined time-dependent incident and transmission monitors. Note that the last run is the depolarized run:

mons = sans_workflow.compute(
    (
        RawMonitor[TransmissionRun[SampleRun], Incident],
        RawMonitor[TransmissionRun[SampleRun], Transmission],
    )
)
mons = sc.DataGroup(
    incident=mons[RawMonitor[TransmissionRun[SampleRun], Incident]],
    transmission=mons[RawMonitor[TransmissionRun[SampleRun], Transmission]],
)
display(sc.DataGroup(sc.collapse(mons['incident'], keep='tof')).plot())
display(sc.DataGroup(sc.collapse(mons['transmission'], keep='tof')).plot())

The task graph for computing the transmission fraction is:

sans_workflow.visualize(TransmissionFraction[SampleRun], graph_attr={'rankdir': 'LR'})

Compute transmission fractions

There are multiple files which together define the time-dependence of the analyzer cell transmission. Note that as before the final run (time) is the depolarized run:

raw_transmission = sans_workflow.compute(TransmissionFraction[SampleRun])

We can plot the computed transmission fractions:

transmission_depolarized = raw_transmission['time', -1].copy()
transmission = raw_transmission['time', :-1].copy()
trans = sc.DataGroup(
    {f"{time:c}": transmission['time', time] for time in transmission.coords['time']}
)
trans[f'depolarized'] = transmission_depolarized
display(trans.plot())

# Sanity check: Where can cosh yield values that can be fitted?
transmission_empty_glass = 0.9 * sc.Unit('dimensionless')
wavelength = sc.midpoints(transmission.coords['wavelength'])
opacity0 = 0.8797823016804095 * sc.Unit('1/angstrom')
(
    sc.acosh(
        sc.values(transmission)
        * sc.exp(opacity0 * wavelength)
        / transmission_empty_glass
    )
    / (opacity0 * wavelength)
).plot()

# TODO Which wavelength bounds should be used?
wav_min = 2.2 * sc.Unit('angstrom')
wav_max = 2.8 * sc.Unit('angstrom')
transmission_truncated = raw_transmission['wavelength', wav_min:wav_max]
transmission_depolarized = transmission_truncated['time', -1].copy()
transmission = transmission_truncated['time', :-1].copy()

We can now setup the polarization analysis workflow. The previously computed transmission fractions are used as workflow inputs:

he3_workflow = pol.he3.He3CellWorkflow(in_situ=False, incoming_polarized=True)
# TODO Is plus correct here, this is period 0? Do we also have minus data?
he3_workflow[pol.he3.He3AnalyzerTransmissionFractionParallel] = transmission
# TODO Fake empty transmission for now, would need to load different period
he3_workflow[pol.he3.He3AnalyzerTransmissionFractionAntiParallel] = transmission[
    'time', 0:0
]
he3_workflow[
    pol.he3.He3CellTransmissionFractionIncomingUnpolarized[
        pol.Analyzer, pol.Depolarized
    ]
] = transmission_depolarized

# When in_situ=False, these params are used as starting guess for the fit
he3_workflow[pol.he3.He3CellLength[pol.Analyzer]] = 0.1 * sc.Unit('m')
he3_workflow[pol.he3.He3CellPressure[pol.Analyzer]] = 1.0 * sc.Unit('bar')
he3_workflow[pol.he3.He3CellTemperature[pol.Analyzer]] = 300.0 * sc.Unit('K')

he3_workflow[pol.he3.He3TransmissionEmptyGlass[pol.Analyzer]] = transmission_empty_glass
he3_workflow.visualize(
    pol.TransmissionFunction[pol.Analyzer], graph_attr={'rankdir': 'LR'}
)

The workflow can compute the transmission function:

func = he3_workflow.compute(pol.TransmissionFunction[pol.Analyzer])

We can evaluate this transmission function at desired time and wavelength points:

wavelength = sc.linspace('wavelength', start=2, stop=16.0, num=141, unit='angstrom')
time = sc.linspace('time', start=0, stop=100000, num=101, unit='s')
display(func.opacity_function(wavelength=wavelength).plot())
display(func.polarization_function(time=time).plot())
display(func(wavelength=wavelength, time=time, plus_minus='plus').plot(norm='log'))
display(func(wavelength=wavelength, time=time, plus_minus='minus').plot(norm='log'))

trans = func(wavelength=wavelength, time=time, plus_minus='plus')
sc.DataGroup(
    {f"{time:c}": trans['time', time] for time in trans.coords['time'][::20]}
).plot(norm='linear', linestyle='solid', marker=None)

trans = func(wavelength=wavelength, time=time, plus_minus='plus')
sc.DataGroup(
    {f"{wav:c}": trans['wavelength', wav] for wav in trans.coords['wavelength'][::20]}
).plot(norm='log', linestyle='solid', marker=None)

func.opacity_function.opacity0

func.polarization_function.C

func.polarization_function.T1

Correction workflow

In the previous section we have setup the workflow for the analyzer. We also computed the transmission function there, but in production this will be done implicitly by running the entire workflow we will setup here. We can combine this with the workflow for the polarizer to obtain the full correction workflow:

supermirror_workflow = pol.SupermirrorWorkflow()
supermirror_workflow.visualize(pol.TransmissionFunction[pol.Polarizer])

We will use a second-order polynomial supermirror efficiency function:

# Note that these coefficients are meaningless, please fill in correct values!
supermirror_workflow[pol.SupermirrorEfficiencyFunction[pol.Polarizer]] = (
    pol.SecondDegreePolynomialEfficiency(
        a=0.5 * sc.Unit('1/angstrom**2'),
        b=0.4 * sc.Unit('1/angstrom'),
        c=0.3 * sc.Unit('dimensionless'),
    )
)

workflow = pol.PolarizationAnalysisWorkflow(
    polarizer_workflow=supermirror_workflow,
    analyzer_workflow=he3_workflow,
)

For a single channel, the complete workflow looks as follows:

workflow.visualize(
    pol.PolarizationCorrectedData[pol.Up, pol.Up], graph_attr={'rankdir': 'LR'}
)