Absorption correction example#

Assuming single scattering, the neutron intensity per unit solid angle at the position \(\mathbf{p}\) is modeled as

\[I(\lambda, \mathbf{x}) = \int_{sample} S(\lambda, \widehat{\mathbf{p}-\mathbf{x}}) T(\lambda, \mathbf{p}, \mathbf{x}) \ d\mathbf{x}\]

where \(S(\lambda, \hat{\mathbf{s}})\) is the probability that a neutron with wavelength \(\lambda\) scatters in the direction \(\hat{\mathbf{s}}\) and \(T(\lambda, \mathbf{p}, \mathbf{x})\) is the transmission rate for neutrons scattering at \(\mathbf{x}\) towards \(\mathbf{p}\). Here the \(\widehat{\quad \cdot\quad}\) means that the vector is normalized to length 1.

If the detector is far away from the sample relative to the size of the sample then \(\widehat{\mathbf{p}-\mathbf{x}}\) is close to constant in the integral, or if \(S\) represents inhomogeneous scattering, the \(S\) term does not depend on \(\mathbf{x}\) and can be moved outside of the integtral:

\[I(\lambda, \mathbf{x}) \approx S(\lambda, \widehat{\mathbf{p} -\bar{\mathbf{x}}}) \int_{sample} T(\lambda, \mathbf{p}, \mathbf{x}) \ d\mathbf{x}\]

where \(\bar{\mathbf{x}}\) denotes the center of the sample.

To compute the scattering probabiltiy \(S\) from the intensity \(I\) we need to estimate the term

\[C(\lambda, \mathbf{p}) = \int_{sample} T(\lambda, \mathbf{p}, \mathbf{x}) \ d\mathbf{x}.\]

This is the “absorption correction” term.

The transmission fraction is a function of path length \(L\) of the neutron going through the sample

\[T(\lambda, \mathbf{p}, \mathbf{x}) = \exp{\big(-\mu(\lambda) L(\hat{\mathbf{b}}, \mathbf{p}, \mathbf{x})\big)}\]

where \(\mu\) is material dependent and \(\hat{\mathbf{b}}\) is the direction of the incoming neutron.

The path length through the sample depends on the shape of the sample, the scattering point \(\mathbf{x}\) and the detection position \(\mathbf{p}\).

To compute \(C(\lambda, \mathbf{p})\) you can use

scippneutron.absorption.base.compute_transmission_map(
    shape,
    material,
    beam_direction,
    wavelength,
    detector_position
)

where shape and material are sample properties:

shape defines L
material defines \(\mu\)
beam_direction is \(\hat{\mathbf{b}}\)
wavelength is \(\lambda\)
detector_position is \(\mathbf{p}\)

[1]:

import scipp as sc
import numpy as np

from scippneutron.absorption import compute_transmission_map, Cylinder, Material
from scippneutron.atoms import ScatteringParams


# Create a material with scattering parameters to represent a Vanadium sample
sample_material = Material(scattering_params=ScatteringParams.for_isotope('V'), effective_sample_number_density=sc.scalar(1, unit='1/angstrom**3'))

# Create a sample shape
sample_shape = Cylinder(
    symmetry_line=sc.vector([0, 1, 0]),
    # center_of_base is placed slightly below the origin, so that the center of the cylinder is at the origin
    center_of_base=sc.vector([0, -.5, 0], unit='cm'),
    radius=sc.scalar(1., unit='cm'),
    height=sc.scalar(0.3, unit='cm')
)

# Create detector positions, in this case the detector positions are placed in a sphere around the sample
theta = sc.linspace('theta', 0, np.pi, 60, endpoint=True, unit='rad')
phi = sc.linspace('phi', 0, 2 * np.pi, 120, endpoint=False, unit='rad')
r = sc.array(dims='r', values=[5.], unit='m')
dims = (*r.dims, *theta.dims, *phi.dims)

# Detector positions in grid on a sphere around the origin
detector_positions = sc.vectors(
    dims=dims,
    values=sc.concat(
        [
            r * sc.sin(theta) * sc.cos(phi),
            sc.broadcast(r * sc.cos(theta), sizes={**r.sizes, **theta.sizes, **phi.sizes}),
            r * sc.sin(theta) * sc.sin(phi),
        ],
        dim='row',
    )
    .transpose([*dims, 'row'])
    .values,
    unit=r.unit,
)

def transmission_fraction(quadrature_kind):
    'Evaluate the transmission fraction using the selected quadrature kind'
    da = compute_transmission_map(
        sample_shape,
        sample_material,
        beam_direction=sc.vector([0, 0, 1]),
        wavelength=sc.linspace('wavelength', 0.5, 2.5, 10, unit='angstrom'),
        detector_position=detector_positions,
        quadrature_kind=quadrature_kind,
    )
    da.coords['phi'] = phi
    da.coords['theta'] = theta
    return da


def show_correction_map(da):
    'Plotting utility'
    return (
        da['wavelength', 0]['r', 0].plot() /
        da['wavelength', 0]['r', 0]['theta', da.sizes['theta']//2].plot() /
        da['wavelength', 0]['r', 0]['theta', da.sizes['theta']//2]['phi', da.sizes['phi']//4 - da.sizes['phi']//6:da.sizes['phi']//4 + da.sizes['phi']//6].plot() /
        da['wavelength', 0]['r', 0]['theta', da.sizes['theta']//2]['phi', da.sizes['phi']//2 - da.sizes['phi']//6:da.sizes['phi']//2 + da.sizes['phi']//6].plot() /
        da['wavelength', 0]['r', 0]['phi', da.sizes['phi']//2].plot()
    )

[2]:

transmission_fraction('medium')

[3]:

show_correction_map(transmission_fraction('cheap'))

[3]:

../_images/user-guide_absorption-correction_3_0.svg

[4]:

show_correction_map(transmission_fraction('medium'))

[4]:

../_images/user-guide_absorption-correction_4_0.svg

[5]:

show_correction_map(transmission_fraction('expensive'))

[5]:

../_images/user-guide_absorption-correction_5_0.svg

What quadrature should I use?#

Use medium first, if it’d not good enough try expensive.