From 68f671f447796aa6fb1e85948e23ffb65030b8f8 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Thu, 2 Oct 2025 18:26:05 +0200 Subject: [PATCH] Perform reduction action for normalization with less steps directly in reducer class, prepare same for datasets --- libeos/normalisation.py | 52 +++++++++++++++--------------- libeos/reduction.py | 70 ++++++++++++++++++++++++++++++----------- 2 files changed, 77 insertions(+), 45 deletions(-) diff --git a/libeos/normalisation.py b/libeos/normalisation.py index 23818ea..8c8a6d6 100644 --- a/libeos/normalisation.py +++ b/libeos/normalisation.py @@ -13,25 +13,25 @@ from .options import NormalisationMethod from .instrument import Grid -class Normalisation: - normFileList = List[str] - normAngle: float - normMonitor: float - norm_lz: np.ndarray +class LZNormalisation: + file_list = List[str] + angle: float + monitor: float + norm: np.ndarray def __init__(self, reference:EventDatasetProtocol, normalisationMethod: NormalisationMethod, grid: Grid): - self.normAngle = reference.geometry.nu-reference.geometry.mu + self.angle = reference.geometry.nu-reference.geometry.mu lamda_e = reference.data.events.lamda detZ_e = reference.data.events.detZ - self.normMonitor = np.sum(reference.data.pulses.monitor) + 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_lz = np.flip(norm_lz, 1) + self.norm = np.flip(norm_lz, 1) else: # correct for reference sm reflectivity lamda_l = grid.lamda() - theta_z = self.normAngle+reference.geometry.delta_z + theta_z = self.angle+reference.geometry.delta_z lamda_lz = (grid.lz().T*lamda_l[:-1]).T theta_lz = grid.lz()*theta_z qz_lz = 4.0*np.pi*np.sin(np.deg2rad(theta_lz))/lamda_lz @@ -39,36 +39,36 @@ class Normalisation: Rsm_lz = np.ones(np.shape(qz_lz)) Rsm_lz = np.where(qz_lz>0.0217, 1-(qz_lz-0.0217)*(0.0625/0.0217), Rsm_lz) Rsm_lz = np.where(qz_lz>0.0217*5, np.nan, Rsm_lz) - self.norm_lz = norm_lz/Rsm_lz - self.normFileList = [os.path.basename(entry) for entry in reference.file_list] + self.norm = norm_lz/Rsm_lz + self.file_list = [os.path.basename(entry) for entry in reference.file_list] @classmethod - def from_file(cls, filename) -> 'Normalisation': + def from_file(cls, filename) -> 'LZNormalisation': logging.warning(f'normalisation matrix: found and using {filename}') self = super().__new__(cls) with open(filename, 'rb') as fh: - self.normFileList = np.load(fh, allow_pickle=True) - self.normAngle = np.load(fh, allow_pickle=True) - self.norm_lz = np.load(fh, allow_pickle=True) - self.normMonitor = np.load(fh, allow_pickle=True) + self.file_list = np.load(fh, allow_pickle=True) + self.angle = np.load(fh, allow_pickle=True) + self.norm = np.load(fh, allow_pickle=True) + self.monitor = np.load(fh, allow_pickle=True) return self @classmethod - def unity(cls, grid:Grid) -> 'Normalisation': + def unity(cls, grid:Grid) -> 'LZNormalisation': logging.warning(f'normalisation is unity') self = super().__new__(cls) - self.norm_lz = grid.lz() - self.normFileList = [] - self.normAngle = 1. - self.normMonitor = 1. + self.norm = grid.lz() + self.file_list = [] + self.angle = 1. + self.monitor = 1. return self def safe(self, filename): with open(filename, 'wb') as fh: - np.save(fh, np.array(self.normFileList), allow_pickle=False) - np.save(fh, np.array(self.normAngle), allow_pickle=False) - np.save(fh, self.norm_lz, allow_pickle=False) - np.save(fh, self.normMonitor, allow_pickle=False) + np.save(fh, np.array(self.file_list), allow_pickle=False) + np.save(fh, np.array(self.angle), allow_pickle=False) + np.save(fh, self.norm, allow_pickle=False) + np.save(fh, self.monitor, allow_pickle=False) def update_header(self, header:Header): - header.measurement_additional_files = self.normFileList + header.measurement_additional_files = self.file_list diff --git a/libeos/reduction.py b/libeos/reduction.py index c189e20..6f02305 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -5,12 +5,13 @@ import sys import numpy as np from orsopy import fileio -from .file_reader import AmorData +from .file_reader import AmorData, AmorEventData from .header import Header from .path_handling import PathResolver from .options import EOSConfig, IncidentAngle, MonitorType, NormalisationMethod from .instrument import Grid -from .normalisation import Normalisation +from .normalisation import LZNormalisation +from . import event_handling as eh, event_analysis as ea MONITOR_UNITS = { MonitorType.neutron_monitor: 'cnts', @@ -42,6 +43,38 @@ class AmorReduction: self.path_resolver = PathResolver(self.reader_config.year, self.reader_config.rawPath) + # setup all actions performed on event datasets before projection on the grid + # The order of these corrections matter as some rely on parameters modified before + if self.reduction_config.normalisationFileIdentifier: + # explicit steps performed on AmorEventDataset for normalization files + self.normevent_actions = eh.ApplyPhaseOffset(self.experiment_config.chopperPhaseOffset) + self.normevent_actions |= eh.CorrectChopperPhase() + self.normevent_actions |= eh.AssociatePulseWithMonitor(self.experiment_config.monitorType, + self.experiment_config.lowCurrentThreshold) + self.normevent_actions |= eh.FilterStrangeTimes() + self.normevent_actions |= eh.MergeFrames() + self.normevent_actions |= ea.AnalyzePixelIDs(self.experiment_config.yRange) + self.normevent_actions |= ea.TofTimeCorrection(self.experiment_config.incidentAngle==IncidentAngle.alphaF) + self.normevent_actions |= ea.CalculateWavelength(self.experiment_config.lambdaRange) + self.normevent_actions |= ea.ApplyMask() + # Actions on datasets not used for normalization + self.dataevent_actions = eh.ApplyPhaseOffset(self.experiment_config.chopperPhaseOffset) + self.dataevent_actions |= eh.ApplyParameterOverwrites(self.experiment_config) # some actions use instrument parameters, change before that + self.dataevent_actions |= eh.CorrectChopperPhase() + self.dataevent_actions |= eh.ExtractWalltime() + self.dataevent_time_correction = eh.CorrectSeriesTime(0) # will be set from first dataset + self.dataevent_actions |= self.dataevent_time_correction + self.dataevent_actions |= eh.AssociatePulseWithMonitor(self.experiment_config.monitorType, + self.experiment_config.lowCurrentThreshold) + self.dataevent_actions |= eh.FilterStrangeTimes() + self.dataevent_actions |= eh.MergeFrames() + self.dataevent_actions |= ea.AnalyzePixelIDs(self.experiment_config.yRange) + self.dataevent_actions |= ea.TofTimeCorrection(self.experiment_config.incidentAngle==IncidentAngle.alphaF) + self.dataevent_actions |= ea.CalculateWavelength(self.experiment_config.lambdaRange) + self.dataevent_actions |= ea.CalculateQ(self.experiment_config.incidentAngle) + self.dataevent_actions |= ea.FilterQzRange(self.experiment_config.qzRange) + self.dataevent_actions |= ea.ApplyMask() + self.grid = Grid(self.reduction_config.qResolution, self.reduction_config.qzRange) def reduce(self): @@ -53,7 +86,7 @@ class AmorReduction: if self.reduction_config.normalisationFileIdentifier: self.create_normalisation_map(self.reduction_config.normalisationFileIdentifier[0]) else: - self.norm = Normalisation.unity(self.grid) + self.norm = LZNormalisation.unity(self.grid) # load R(q_z) curve to be subtracted: if self.reduction_config.subtract: @@ -97,7 +130,7 @@ class AmorReduction: self.monitor = np.sum(self.file_reader.monitorPerPulse) logging.warning(f' monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.experiment_config.monitorType]}') qz_lz, qx_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, self.mask_lz = self.project_on_lz( - self.file_reader, self.norm.norm_lz, self.norm.normAngle, lamda_e, detZ_e) + self.file_reader, self.norm.norm, self.norm.angle, lamda_e, detZ_e) #if self.monitor>1 : # ref_lz /= self.monitor # err_lz /= self.monitor @@ -117,7 +150,7 @@ class AmorReduction: # projection on qz-grid q_q, R_q, dR_q, dq_q = self.project_on_qz(qz_lz, ref_lz, err_lz, res_lz, - self.norm.norm_lz, self.mask_lz) + self.norm.norm, self.mask_lz) # The filtering is now done by restricting the qz-grid #filter_q = np.where((self.experiment_config.qzRange[0]>q_q) & (q_q>self.experiment_config.qzRange[1]), @@ -178,7 +211,7 @@ class AmorReduction: lindex_lz.T, tindex_lz.T, int_lz.T, - self.norm.norm_lz.T, + self.norm.norm.T, np.where(self.mask_lz, 1, 0).T, qx_lz.T, ): @@ -215,14 +248,14 @@ class AmorReduction: logging.info(f' {ti:<4d} {time:6.0f} {self.monitor:7.2f} {MONITOR_UNITS[self.experiment_config.monitorType]}') qz_lz, qx_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz = self.project_on_lz( - self.file_reader, self.norm.norm_lz, self.norm.normAngle, lamda_e, detZ_e) + self.file_reader, self.norm.norm, self.norm.angle, lamda_e, detZ_e) try: ref_lz *= self.reduction_config.scale[i] err_lz *= self.reduction_config.scale[i] except IndexError: ref_lz *= self.reduction_config.scale[-1] err_lz *= self.reduction_config.scale[-1] - q_q, R_q, dR_q, dq_q = self.project_on_qz(qz_lz, ref_lz, err_lz, res_lz, self.norm.norm_lz, mask_lz) + q_q, R_q, dR_q, dq_q = self.project_on_qz(qz_lz, ref_lz, err_lz, res_lz, self.norm.norm, mask_lz) filter_q = np.where((self.experiment_config.qzRange[0] 1e6: self.norm.safe(n_path) @@ -444,7 +476,7 @@ class AmorReduction: logging.error('unknown normalisation method! Use [u]nder, [o]ver or [d]irect illumination') ref_lz = (int_lz / norm_lz) if self.monitor > 1e-6 : - ref_lz *= self.norm.normMonitor / self.monitor + ref_lz *= self.norm.monitor/self.monitor else: logging.info(' low monitor -> nan output') ref_lz *= np.nan