stitch - Overlapping reflectivity curve stitching

poisson_average

Compute the poisson average of y/dy using a set of data points.

stitch

Stitch together multiple measurements into one.

Data stitching.

Join together datasets yielding unique sorted x.

refl1d.stitch.poisson_average(xdxydyw)[source]

Compute the poisson average of y/dy using a set of data points.

The returned x, dx is the weighted average of the inputs:

x = sum(x*I)/sum(I)
dx = sum(dx*I)/sum(I)

The returned y, dy use Poisson averaging:

w = sum(y/dy^2)
y = sum((y/dy)^2)/w
dy = sqrt(y/w)

The above formula gives the expected result for combining two measurements, assuming there is no uncertainty in the monitor.

measure N counts during M monitors
  rate:                   r = N/M
  rate uncertainty:       dr = sqrt(N)/M
  weighted rate:          r/dr^2 = (N/M) / (N/M^2) =  M
  weighted rate squared:  r^2/dr^2 = (N^2/M^2) / (N/M^2) = N

for two measurements Na, Nb
  w = ra/dra^2 + rb/drb^2 = Ma + Mb
  y = ((ra/dra)^2 + (rb/drb)^2)/w = (Na + Nb)/(Ma + Mb)
  dy = sqrt(y/w) = sqrt( (Na + Nb)/ w^2 ) = sqrt(Na+Nb)/(Ma + Mb)

This formula isn’t strictly correct when applied to values which have been scaled, for example to account for an attenuator in the counting system.

refl1d.stitch.stitch(data, same_x=0.001, same_dx=0.001)[source]

Stitch together multiple measurements into one.

data a list of datasets with x, dx, y, dy attributes same_x minimum point separation (default is 0.001). same_dx minimum change in resolution that may be averaged (default is 0.001).

WARNING: the returned x values may be data dependent, with two measured sets having different x after stitching, even though the measurement conditions are identical!!

Either add an intensity weight to the datasets:

probe.I = slitscan

or use interpolation if you need to align two stitched scans:

import numpy as np
x1, dx1, y1, dy1 = stitch([a1, b1, c1, d1])
x2, dx2, y2, dy2 = stitch([a2, b2, c2, d2])
x2[0], x2[-1] = x1[0], x1[-1] # Force matching end points
y2 = np.interp(x1, x2, y2)
dy2 = np.interp(x1, x2, dy2)
x2 = x1

WARNING: the returned dx value underestimates the true x, depending on the relative weights of the averaged data points.