Possible error in reflectometry footprint correction

[jokasimr]

The Amor reduction uses a "footprint" correction to scale the number of counts on the detector. From what I understand, the footprint correction compensates for the size of the beam on sample depending on the incidence angle $\theta$ and perhaps extending past the sample, reducing the total number of neutrons hitting the sample from small angles.

The correction is defined in src/ess/reflectometry/corrections.py. I think what it does is it computes the fraction of the beam intensity hitting the sample based on the assumptions that

1. the intensity of the beam is well approximated by a normal distribution over the cross section,
2. the beam hits the center of the sample,
3. the sample width is much larger than the FWHM of the beam.

Let the full width of the sample be $w$, the FWHM of the beam be $d$ and the FWHM of the standard normal distribution be $d_0$. From what I can tell the correction factor used in the code is $\mathrm{erf} \big( \frac{d_0 w \sin(\theta)}{d}\big)$, if sample_size in the code corresponds to $w$. But I think the correction should be $\mathrm{erf} \big( \frac{d_0 w \sin(\theta)}{2\sqrt{2}d}\big)$.

Why? A neutron in the beam at a distance $r$ from the center hits the sample if $\frac{r}{\sin{\theta}} < \frac{w}{2}$. The fraction of the intensity that hits the sample is

$I{hit} = \int{-\infty}^{\infty} \mathcal{1}_{\frac{r}{\sin{\theta}} < \frac{w}{2}} I(r) \ dr $

$= \int_{-\sin{\theta} w / 2}^{\sin{\theta} w / 2} I(r) \ dr$

$= 2 I{total} \int{0}^{\sin{\theta} w / 2} N(r; 0, d/d_0) \ dr $

$= 2 I{total} \int{0}^{\frac{d_0 \sin{\theta} w}{2 d}} N(x ; 0, 1) \ dx $

$= I_{total} \ \mathrm{erf}\bigg({\frac{w}{2}\frac{d_0 \sin{\theta}}{\sqrt{2} d}}\bigg)$.

Here $N(x; \mu, \sigma)$ denotes a normal distribution with standard deviation $\sigma$.

(see Python script for verifying the integration)

import sympy as sp
r, d, d0 = sp.symbols('r, d, d_0', real=True)
theta, sigma, w = sp.symbols('\\theta, \sigma, w', positive=True)

# Bell curve with standard deviation sigma
G = sp.exp(-r**2 / (2 * sigma**2)) / (sigma * sp.sqrt(2 * sp.pi))
# Verify integral is 1
assert sp.integrate(G, (r, '-oo', 'oo')) == 1

inside_of_sample = sp.Piecewise(
    (0, 2 * sp.Abs(r) > sp.Abs(sp.sin(theta)) * w),
    (1, True),
)
sp.integrate(G * inside_of_sample, (r, '-oo', 'oo')).subs({sigma: d/d0})

Part of the discrepancy might be that sample_size actually denotes the radius of the sample or "half the width" of the sample if it is square. What should it be? However, there's still a factor $\sqrt{2}$ unaccounted for.

Might just be I've misunderstood something.

Didn't find the description of this correction procedure in the reference in the workflow.

[SimonHeybrock]

@arm61 Do you know if there is a reference where this equation was taken from?

[SimonHeybrock]

@jokasimr What about this factor: https://github.com/scipp/ess/blob/7025e1770239899394a75e33fe3dfa997ff92649/src/ess/amor/tools.py#L8-L11 ... looks similar to what you have, so the code I see does not look like the equation you wrote?

[jokasimr]

_STD_TO_FWHM should be equal to $d_0$ in the description above.

[SimonHeybrock]

Did you assume a spherical or a square/rectangular beam profile? Does it make a difference?

[jokasimr]

I assumed the intensity over the cross section of the beam is described by $I(r, t) = I_{total} N(r; 0, d/d_0) \ \frac{S(t/d)}{d}$ where $S(x)=1$ ifor $x\in[-1/2, 1/2]$ and zero otherwise (square profile). $\hat{\mathbf{t}}$ is a direction that is parallel to both the beam cross section and to the sample plane. $\hat{\mathbf{r}}$ is parallel to the beam cross section but orthogonal to $\hat{\mathbf{t}}$.

I think as long as $d < < w$ it doesn't matter much what shape we prescribe in the direction $\hat{\mathbf{t}}$, so to simplify I removed it and treated the beam cross section as 1-dimensional.

[SimonHeybrock]

I think as long as d<<w it doesn't matter much what shape we prescribe in the direction t^, so to simplify I removed it and treated the beam cross section as 1-dimensional.

Not so sure about this: Say the beam profile is roughly circular (with Gaussian drop off), so it illuminates an "elliptical" area of the sample plane. If you consider the center of the beam, then your math applies. But if we move to the left or right (your $\hat t $, I think?), the distribution is "narrower" so the fraction that falls outside the sample is smaller?

[arm61]

This is the paper that most of the logic in that reduction come from https://doi.org/10.1016/j.nima.2016.03.007

[jokasimr]

But if we move to the left or right (your t, I think?), the distribution is "narrower" so the fraction that falls outside the sample is smaller?

Good point. I think the 1D approach used above to find the correction factor is only right if either

• the intensity is constant in $\hat{\mathbf{t}}$ and $d < w$,
• or if $d/w \to 0$ so that the distribution in $\hat{\mathbf{t}}$ is narrow like a delta pulse.

[jokasimr]

It seems quite unclear how this should actually be done for Amor. For ESTIA it might not even be needed. So for now I'm removing this from selected.