Perform reduction action for normalization with less steps directly in reducer class, prepare same for datasets

This commit is contained in:
2025-10-02 18:26:05 +02:00
parent cb4415ad3d
commit 68f671f447
2 changed files with 77 additions and 45 deletions
+26 -26
View File
@@ -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
+51 -19
View File
@@ -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]<q_q) & (q_q<self.experiment_config.qzRange[1]),
True, False)
@@ -344,22 +377,21 @@ class AmorReduction:
def create_normalisation_map(self, short_notation):
outputPath = self.output_config.outputPath
normalisation_list = self.path_resolver.expand_file_list(short_notation)
name = str(normalisation_list[0])
for norm_i in normalisation_list[1:]:
name += f'_{norm_i}'
name = '_'.join(map(str, normalisation_list))
n_path = os.path.join(outputPath, f'{name}_{str(self.experiment_config.monitorType)}.norm')
if os.path.exists(n_path):
logging.warning(f'normalisation matrix: found and using {n_path}')
self.norm = Normalisation.from_file(n_path)
self.header.measurement_additional_files = self.norm.normFileList
self.norm = LZNormalisation.from_file(n_path)
self.header.measurement_additional_files = self.norm.file_list
else:
logging.warning(f'normalisation matrix: using the files {normalisation_list}')
reference = AmorData(header=self.header,
reader_config=self.reader_config,
config=self.experiment_config,
short_notation=short_notation, norm=True).dataset
self.norm = Normalisation(reference, self.reduction_config.normalisationMethod, self.grid)
normalization_files = list(map(self.path_resolver.get_path, normalisation_list))
reference = AmorEventData(normalization_files[0])
for nfi in normalization_files[1:]:
reference.append(AmorEventData(nfi))
self.normevent_actions(reference)
self.norm = LZNormalisation(reference, self.reduction_config.normalisationMethod, self.grid)
if reference.data.events.shape[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