sinq-reflectometry/eos
0
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 5 Wiki Activity

Re-write qz projection to work with 0-count norm bins, don't filter norm by counts as default, optional smoothing of norm, fix dQ calculation

Browse Source
This commit is contained in:
glavic_a
2025-12-01 15:11:49 +01:00
parent 6597ca22d1
commit a56f29191d
4 changed files with 75 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files
eos/normalization.py
+29 -1
View File
@@ -24,7 +24,6 @@ class LZNormalisation:
detZ_e = reference.data.events.detZ
self.monitor = np.sum(reference.data.pulses.monitor)
norm_lz, _, _ = np.histogram2d(lamda_e, detZ_e, bins=(grid.lamda(), grid.z()))
norm_lz = np.where(norm_lz>2, norm_lz, np.nan)
if normalisationMethod==NormalisationMethod.direct_beam:
self.norm = np.flip(norm_lz, 1)
else:
@@ -100,3 +99,32 @@ class LZNormalisation:
def update_header(self, header:Header):
header.measurement_additional_files = self.file_list
def smooth(self, width):
logging.debug(f'apply convolution with gaussian along lambda axis to smooth norm, sigma={width}')
try:
from scipy.signal import fftconvolve
except ImportError:
self._smooth_slow(width)
kx = np.arange(self.norm.shape[0])-self.norm.shape[0]/2.
kernel = np.exp(-0.5*kx**2/width**2)
kernel/=kernel.sum()
kernel = kernel[:, np.newaxis]*np.ones(self.norm.shape[1])[np.newaxis, :]
unorm = np.where(np.isnan(self.norm), 0., self.norm)
nnorm = fftconvolve(unorm, kernel, mode='same', axes=0)
nnorm[np.isnan(self.norm)] = np.nan
self.norm = nnorm
def _smooth_slow(self, width):
# like smooth but using numpy buildin slow convolve
nnorm = np.zeros_like(self.norm)
kx = np.arange(self.norm.shape[0])-self.norm.shape[0]/2.
kernel = np.exp(-0.5*kx**2/width**2)
kernel/=kernel.sum()
unorm = np.where(np.isnan(self.norm), 0., self.norm)
for row in range(self.norm.shape[1]):
nnorm[:, row] = np.convolve(unorm[:, row], kernel, mode='same')
nnorm[np.isnan(self.norm)] = np.nan
self.norm = nnorm
eos/options.py
+15
View File
@@ -381,6 +381,21 @@ class ReflectivityReductionConfig(ArgParsable):
'priority': 90,
'group': 'input data',
'help': 'normalisation method: [o]verillumination, [u]nderillumination, [d]irect_beam'})
normalizationFilter: float = field(
default=-1,
metadata={
'group': 'input data',
'help': 'minimum normalization counts in lambda-theta bin to use, else filter'})
normAngleFilter: float = field(
default=-1,
metadata={
'group': 'input data',
'help': 'minimum normalization counts total thetat bin to use, else filter'})
normalizationSmoothing: float = field(
default=0,
metadata={
'group': 'input data',
'help': 'apply convolution on lambda axes to smooth the normalization data, useful for low statistics'})
scale: List[float] = field(
default_factory=lambda: [1.],
metadata={
eos/projection.py
+26 -24
View File
@@ -176,9 +176,12 @@ class LZProjection(ProjectionInterface):
self.data.mask &= self.lamda>=lamda_range[0]
self.data.mask &= self.lamda<=lamda_range[1]
def apply_norm_mask(self, norm: LZNormalisation):
def apply_norm_mask(self, norm: LZNormalisation, min_norm=-1, min_theta=-1):
# Mask points where normliazation is nan
self.data.mask &= np.logical_not(np.isnan(norm.norm))
self.data.mask &= np.logical_not(np.isnan(norm.norm))&(norm.norm>min_norm)
if min_theta>0:
thsum = np.nansum(norm.norm, axis=0)
self.data.mask &= (thsum>min_theta)[np.newaxis, :]
def project(self, dataset: EventDatasetProtocol, monitor: float):
"""
@@ -198,6 +201,7 @@ class LZProjection(ProjectionInterface):
self.data[:] = 0
self.data.mask = True
self.monitor = 0.
self.norm_monitor = 1.
@property
def I(self):
@@ -207,7 +211,7 @@ class LZProjection(ProjectionInterface):
def calc_error(self):
# calculate error bars for resulting intensity after normalization
self.data.err = self.data.ref * np.sqrt( 1./(self.data.I+.1) + 1./self.data.norm )
self.data.err = self.data.ref * np.sqrt( 1./(self.data.I+.1) + 1./(self.data.norm+0.1) )
def normalize_over_illuminated(self, norm: LZNormalisation):
"""
@@ -221,10 +225,11 @@ class LZProjection(ProjectionInterface):
fp_corr_lz = np.where(np.absolute(delta_lz+norm.angle)>5e-3,
(delta_lz+self.angle)/(delta_lz+norm.angle), np.nan)
self.data.mask &= np.logical_not(np.isnan(fp_corr_lz))
ref_lz = self.data.I/norm_lz/fp_corr_lz
self.data.norm = norm_lz*fp_corr_lz
self.norm_monitor = norm.monitor
ref_lz = self.data.I/np.where(self.data.norm>0, self.data.norm, np.nan)
ref_lz *= norm.monitor/self.monitor
ref_lz[np.logical_not(self.data.mask)] = np.nan
self.data.norm = norm_lz
self.data.ref = ref_lz
self.calc_error()
self.is_normalized = True
@@ -244,6 +249,7 @@ class LZProjection(ProjectionInterface):
raise ValueError("Dataset needs to be normalized, first")
self.data.ref *= factor
self.data.err *= factor
self.norm_monitor /= factor
def project_on_qz(self):
if not self.is_normalized:
@@ -251,29 +257,25 @@ class LZProjection(ProjectionInterface):
q_q = self.grid.q()
weights_lzf = self.data.norm[self.data.mask]
q_lzf = self.data.qz[self.data.mask]
R_lzf = self.data.ref[self.data.mask]
dR_lzf = self.data.err[self.data.mask]
I_lzf = self.data.I[self.data.mask]
dq_lzf = self.data.res[self.data.mask]
N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf )[0]
# get number of grid points contributing to a bin, filter points with no contribution
N_q = np.histogram(q_lzf, bins = q_q)[0]
N_q = np.where(N_q > 0, N_q, np.nan)
fltr = N_q>0
R_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf * R_lzf )[0]
R_q = R_q / N_q
dR_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dR_lzf)**2 )[0]
dR_q = np.sqrt( dR_q ) / N_q
# TODO: different error propagations for dR and dq!
# this is what should work:
#dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0]
#dq_q = np.sqrt( dq_q ) / N_q
# and this actually works:
N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf**2 )[0]
N_q = np.where(N_q > 0, N_q, np.nan)
dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0]
dq_q = np.sqrt( dq_q / N_q )
# calculate sum of all normalization weights per bin
W_q = np.maximum(np.histogram(q_lzf, bins = q_q, weights = weights_lzf)[0], 1e-10)
# calculate sum of all dataset counts per bin
I_q = np.histogram(q_lzf, bins = q_q, weights = I_lzf)[0]
# normlaize dataaset by normalization counts and scale by monitor
R_q = np.where(fltr, I_q*self.norm_monitor/self.monitor / W_q, np.nan)
# error as squar-root of counts and sqrt from normalization (dR/R = sqrt( (dI/I)² + (dW/W)²)
dR_q = np.where(fltr, R_q*(np.sqrt(1./(I_q+0.1)+ 1./(W_q+0.1))), np.nan)
# q-resolution is the weighted sum of individual points q-resolutions
dq_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf * dq_lzf**2 )[0]
dq_q = np.where(fltr, np.sqrt(dq_q/W_q), np.nan)
return ProjectedReflectivity(R_q, dR_q, (q_q[1:]+q_q[:-1])/2., dq_q)
def plot(self, **kwargs):
eos/reduction_reflectivity.py
+5 -1
View File
@@ -85,6 +85,9 @@ class ReflectivityReduction:
else:
self.norm = LZNormalisation.unity(self.grid)
if self.config.reduction.normalizationSmoothing:
self.norm.smooth(self.config.reduction.normalizationSmoothing)
# load R(q_z) curve to be subtracted:
if self.config.reduction.subtract:
self.sq_q, self.sR_q, self.sdR_q, self.sFileName = self.loadRqz(self.config.reduction.subtract)
@@ -390,7 +393,8 @@ class ReflectivityReduction:
proj.apply_lamda_mask(self.config.experiment.lambdaRange)
proj.apply_norm_mask(self.norm)
proj.apply_norm_mask(self.norm, min_norm=self.config.reduction.normalizationFilter,
min_theta=self.config.reduction.normAngleFilter)
proj.project(dataset, self.monitor)