From 31201795b0d844c669ee170c68fc32c808c7746e Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 7 Oct 2025 16:19:57 +0200 Subject: [PATCH] implement events2histogram for LZ graph and automatic update --- eos/e2h.py | 12 ++- eos/event_analysis.py | 2 +- eos/event_data_types.py | 2 +- eos/event_handling.py | 3 + eos/file_reader.py | 32 ++++--- eos/options.py | 8 ++ eos/path_handling.py | 34 +++++++- eos/projection.py | 27 +++++- eos/reduction_e2h.py | 153 +++++++++++++++++++++++++++++++--- eos/reduction_reflectivity.py | 2 +- setup.cfg | 1 + 11 files changed, 235 insertions(+), 41 deletions(-) diff --git a/eos/e2h.py b/eos/e2h.py index 5d04981..9d49646 100644 --- a/eos/e2h.py +++ b/eos/e2h.py @@ -1,28 +1,32 @@ import logging # need to do absolute import here as pyinstaller requires it -from eos.options import E2HConfig, ReaderConfig, ExperimentConfig +from eos.options import E2HConfig, ReaderConfig, ExperimentConfig, E2HReductionConfig from eos.command_line import commandLineArgs from eos.logconfig import setup_logging, update_loglevel +from eos.reduction_e2h import E2HReduction + def main(): setup_logging() + logging.getLogger('matplotlib').setLevel(logging.WARNING) # read command line arguments and generate classes holding configuration parameters - clas = commandLineArgs([ReaderConfig, ExperimentConfig], + clas = commandLineArgs([ReaderConfig, ExperimentConfig, E2HReductionConfig], 'events2histogram') update_loglevel(clas.verbose) reader_config = ReaderConfig.from_args(clas) experiment_config = ExperimentConfig.from_args(clas) - config = E2HConfig(reader_config, experiment_config, ) + reduction_config = E2HReductionConfig.from_args(clas) + config = E2HConfig(reader_config, experiment_config, reduction_config) logging.warning('######## events2histogram - data vizualization for Amor ########') # only import heavy module if sufficient command line parameters were provided from eos.reduction_reflectivity import ReflectivityReduction # Create reducer with these arguments - reducer = ReflectivityReduction(config) + reducer = E2HReduction(config) # Perform actual reduction reducer.reduce() diff --git a/eos/event_analysis.py b/eos/event_analysis.py index 250434e..ade1bb6 100644 --- a/eos/event_analysis.py +++ b/eos/event_analysis.py @@ -18,7 +18,7 @@ class ExtractWalltime(EventDataAction): def perform_action(self, dataset: EventDatasetProtocol) ->None: wallTime = extract_walltime(dataset.data.events.tof, dataset.data.packets.start_index, - dataset.data.packets.Time) + dataset.data.packets.time) logging.debug(f' expending event stream by wallTime') new_events = append_fields(dataset.data.events, [('wallTime', wallTime.dtype)]) new_events.wallTime = wallTime diff --git a/eos/event_data_types.py b/eos/event_data_types.py index 2faf297..c4bcb43 100644 --- a/eos/event_data_types.py +++ b/eos/event_data_types.py @@ -31,7 +31,7 @@ class AmorTiming: # Structured datatypes used for event streams EVENT_TYPE = np.dtype([('tof', np.float64), ('pixelID', np.uint32), ('mask', np.int32)]) -PACKET_TYPE = np.dtype([('start_index', np.uint32), ('Time', np.int64)]) +PACKET_TYPE = np.dtype([('start_index', np.uint32), ('time', np.int64)]) PULSE_TYPE = np.dtype([('time', np.int64), ('monitor', np.float32)]) PC_TYPE = np.dtype([('current', np.float32), ('time', np.int64)]) diff --git a/eos/event_handling.py b/eos/event_handling.py index 0b51c6d..ccd39db 100644 --- a/eos/event_handling.py +++ b/eos/event_handling.py @@ -149,6 +149,9 @@ class TofTimeCorrection(EventDataAction): self.correct_chopper_opening = correct_chopper_opening def perform_action(self, dataset: EventDatasetProtocol) ->None: + if not 'delta' in dataset.data.events.dtype.names: + raise ValueError( + "TofTimeCorrection requires delta to be extracted, please run AnalyzePixelIDs first") d = dataset.data if self.correct_chopper_opening: d.events.tof -= ( d.events.delta / 180. ) * dataset.timing.tau diff --git a/eos/file_reader.py b/eos/file_reader.py index a8899e5..5a59c68 100644 --- a/eos/file_reader.py +++ b/eos/file_reader.py @@ -30,7 +30,7 @@ AMOR_LOCAL_TIMEZONE = zoneinfo.ZoneInfo(key='Europe/Zurich') if platform.node().startswith('amor'): NICOS_CACHE_DIR = '/home/amor/nicosdata/amor/cache/' - GREP = '/usr/bin/grep "%s"' + GREP = '/usr/bin/grep "value"' else: NICOS_CACHE_DIR = None @@ -63,10 +63,11 @@ class AmorHeader: try: logging.warning(f" using parameter {nicos_key} from nicos cache") year_date = self.fileDate.strftime('%Y') - value = str(subprocess.getoutput(f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year_date}')).split('\t')[-1] + call = f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year_date}' + value = str(subprocess.getoutput(call)).split('\t')[-1] return dtype(value) except Exception: - logging.error("Couldn't get value from nicos cache", exc_info=True) + logging.error(f"Couldn't get value from nicos cache {nicos_key}, {call}") return dtype(0) else: logging.warning(f" parameter {key} not found, relpace by zero") @@ -151,7 +152,7 @@ class AmorHeader: chopperSeparation, detectorDistance, chopperDetectorDistance) self.timing = AmorTiming(ch1TriggerPhase, ch2TriggerPhase, chopperSpeed, chopperPhase, tau) - polarizationConfigLabel = self._replace_if_missing('polarization/configuration/average_value', 'polarizer_config_label', int) + polarizationConfigLabel = self._replace_if_missing('polarization/configuration/average_value', 'polarization_config_label', int) polarizationConfig = fileio.Polarization(polarizationConfigs[polarizationConfigLabel]) logging.debug(f' polarization configuration: {polarizationConfig} (index {polarizationConfigLabel})') @@ -236,10 +237,6 @@ class AmorEventData(AmorHeader): self.hdf = hdf self.read_event_stream() - # actions applied to any dataset - self.read_chopper_trigger_stream() - self.read_proton_current_stream() - if type(fileName) is str: # close the input file to free memory, only if the file was opened in this object self.hdf.close() @@ -253,12 +250,13 @@ class AmorEventData(AmorHeader): """ packets = np.recarray(self.hdf['/entry1/Amor/detector/data/event_index'].shape, dtype=PACKET_TYPE) packets.start_index = self.hdf['/entry1/Amor/detector/data/event_index'][:] - packets.Time = self.hdf['/entry1/Amor/detector/data/event_time_zero'][:] + packets.time = self.hdf['/entry1/Amor/detector/data/event_time_zero'][:] try: # packet index that matches first event index start_packet = int(np.where(packets.start_index==self.first_index)[0][0]) except IndexError: - raise IndexError(f'No event packet found starting at event #{self.first_index}') + raise EOFError(f'No event packet found starting at event #{self.first_index}, ' + f'number of events is {self.hdf["/entry1/Amor/detector/data/event_time_offset"].shape[0]}') packets = packets[start_packet:] nevts = self.hdf['/entry1/Amor/detector/data/event_time_offset'].shape[0] @@ -279,15 +277,15 @@ class AmorEventData(AmorHeader): events.pixelID = self.hdf['/entry1/Amor/detector/data/event_id'][self.first_index:self.last_index+1] events.mask = 0 - pulses = self.read_chopper_trigger_stream() - current = self.read_proton_current_stream() + pulses = self.read_chopper_trigger_stream(packets) + current = self.read_proton_current_stream(packets) self.data = AmorEventStream(events, packets, pulses, current) if self.first_index>0 and not self.EOF: # label the file name if not all events were used self.file_list[0] += f'[{self.first_index}:{self.last_index}]' - def read_chopper_trigger_stream(self): + def read_chopper_trigger_stream(self, packets): chopper1TriggerTime = np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][:-2], dtype=np.int64) #self.chopper2TriggerTime = self.chopper1TriggerTime + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time'][:-2], dtype=np.int64) # + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][:], dtype=np.int64) @@ -304,17 +302,17 @@ class AmorEventData(AmorHeader): pulses.monitor = 1. # default is monitor pulses as it requires no calculation # apply filter in case the events were filtered if self.first_index>0 or not self.EOF: - pulses = pulses[(pulses.time>=self.data.packets.Time[0])&(pulses.time<=self.data.packets.Time[-1])] + pulses = pulses[(pulses.time>=packets.time[0])&(pulses.time<=packets.time[-1])] self.eventStartTime = startTime return pulses - def read_proton_current_stream(self): + def read_proton_current_stream(self, packets): proton_current = np.recarray(self.hdf['entry1/Amor/detector/proton_current/time'].shape, dtype=PC_TYPE) proton_current.time = self.hdf['entry1/Amor/detector/proton_current/time'][:] proton_current.current = self.hdf['entry1/Amor/detector/proton_current/value'][:,0] if self.first_index>0 or not self.EOF: - proton_current = proton_current[(proton_current.time>=self.data.packets.Time[0])& - (proton_current.time<=self.data.packets.Time[-1])] + proton_current = proton_current[(proton_current.time>=packets.time[0])& + (proton_current.time<=packets.time[-1])] return proton_current def info(self): diff --git a/eos/options.py b/eos/options.py index 2866198..269121d 100644 --- a/eos/options.py +++ b/eos/options.py @@ -661,6 +661,14 @@ class E2HReductionConfig(ArgParsable): }, ) + plot_colormap: PlotColormaps = field( + default=PlotColormaps.gist_ncar, + metadata={ + 'short': 'pcmap', + 'group': 'output', + 'help': 'matplotlib colormap used in lambda-theta graphs when plotting', + }, + ) @dataclass diff --git a/eos/path_handling.py b/eos/path_handling.py index e401d22..2783ddc 100644 --- a/eos/path_handling.py +++ b/eos/path_handling.py @@ -1,6 +1,7 @@ """ Defines how file paths are resolved from short_notation, year and number to filename. """ +import logging import os from typing import List @@ -18,7 +19,7 @@ class PathResolver: """Evaluate string entry for file number lists""" file_list = [] for i in short_notation.split(','): - if '-' in i: + if '-' in i and not i.startswith('-'): if ':' in i: step = i.split(':', 1)[1] file_list += range(int(i.split('-', 1)[0]), @@ -34,6 +35,8 @@ class PathResolver: return list(sorted(file_list)) def get_path(self, number): + if number<=0: + number = self.search_latest(number) fileName = f'amor{self.year}n{number:06d}.hdf' path = '' for rawd in self.rawPath: @@ -41,9 +44,32 @@ class PathResolver: path = rawd break if not path: - if os.path.exists( - f'/afs/psi.ch/project/sinqdata/{self.year}/amor/{int(number/1000)}/{fileName}'): - path = f'/afs/psi.ch/project/sinqdata/{self.year}/amor/{int(number/1000)}' + from glob import glob + potential_file = glob(f'/home/amor/data/{self.year}/*/{fileName}') + if len(potential_file)>0: + path = os.path.dirname(potential_file[0]) else: raise FileNotFoundError(f'# ERROR: the file {fileName} can not be found in {self.rawPath}') return os.path.join(path, fileName) + + def search_latest(self, number): + if number>0: + raise ValueError('number needs to be relative index (negative)') + if os.path.exists(f'/home/amor/data/{self.year}/DataNumber'): + try: + with open(f'/home/amor/data/{self.year}/DataNumber', 'r') as fh: + current_index = int(fh.readline())-1 + except Exception: + logging.error('Can not access DataNumber', exc_info=True) + else: + return current_index+number + # find all files from the given year, convert to number and return latest + from glob import glob + possible_files = [] + for rawd in self.rawPath: + possible_files += glob(os.path.join(rawd, f'amor{self.year}n??????.hdf')) + + possible_files += glob(f'/home/amor/data/{self.year}/*/amor{self.year}n??????.hdf') + possible_indices = list(set([int(os.path.basename(fi)[9:15]) for fi in possible_files])) + possible_indices.sort() + return possible_indices[number-1] diff --git a/eos/projection.py b/eos/projection.py index 894208d..b397d94 100644 --- a/eos/projection.py +++ b/eos/projection.py @@ -3,6 +3,8 @@ Classes used to calculate projections/binnings from event data onto given grids. """ import logging +from typing import Protocol + import numpy as np from dataclasses import dataclass @@ -75,6 +77,12 @@ class ProjectedReflectivity: plt.xlabel('Q / $\\AA^{-1}$') plt.ylabel('R') +class ProjectionInterface(Protocol): + def project(self, dataset: EventDatasetProtocol, monitor: float): ... + def clear(self): ... + def plot(self, **kwargs): ... + def update_plot(self): ... + class LZProjection: grid: LZGrid lamda: np.ndarray @@ -165,6 +173,7 @@ class LZProjection: Project dataset on grid and add to intensity. Can be called multiple times to sequentially add events. """ + # TODO: maybe move monitor calculation in here instead of reduction? e = dataset.data.events int_lz, *_ = np.histogram2d(e.lamda, e.detZ, bins = (self.grid.lamda(), self.grid.z())) self.data.I += int_lz @@ -172,6 +181,11 @@ class LZProjection: # in case the intensity changed one needs to normalize again self.is_normalized = False + def clear(self): + # empty data + self.data[:] = 0 + self.data.mask = True + @property def I(self): output = self.data.I[:] @@ -304,14 +318,23 @@ class LZProjection: if self.is_normalized: if not 'norm' in kwargs: kwargs['norm'] = LogNorm(2e-3, 2.0) - plt.pcolormesh(self.lamda, self.alphaF, self.data.ref, **kwargs) + self._graph = plt.pcolormesh(self.lamda, self.alphaF, self.data.ref, **kwargs) if cmap: plt.colorbar(label='R') else: if not 'norm' in kwargs: kwargs['norm'] = LogNorm() - plt.pcolormesh(self.lamda, self.alphaF, self.data.I, **kwargs) + self._graph = plt.pcolormesh(self.lamda, self.alphaF, self.data.I, **kwargs) if cmap: plt.colorbar(label='I / cpm') plt.xlabel('$\\lambda$ / $\\AA$') plt.ylabel('$\\Theta$ / °') + + def update_plot(self): + """ + Inline update of last plot by just updating the data. + """ + if self.is_normalized: + self._graph.set_array(self.data.ref) + else: + self._graph.set_array(self.data.I) diff --git a/eos/reduction_e2h.py b/eos/reduction_e2h.py index 86203d1..6b93782 100644 --- a/eos/reduction_e2h.py +++ b/eos/reduction_e2h.py @@ -4,19 +4,29 @@ Can be used as a live preview with automatic update when files are modified. """ import logging +import os +import matplotlib.pyplot as plt -from .file_reader import AmorEventData +from orsopy import fileio +from datetime import datetime + +from .file_reader import AmorEventData, AmorHeader from .header import Header from .instrument import LZGrid from .normalization import LZNormalisation -from .options import E2HConfig, IncidentAngle, MonitorType, NormalisationMethod, MONITOR_UNITS, ExperimentConfig +from .options import E2HConfig, E2HPlotArguments, IncidentAngle, MonitorType, E2HPlotSelection from . import event_handling as eh, event_analysis as ea +from .path_handling import PathResolver +from .projection import LZProjection, ProjectionInterface + class E2HReduction: config: E2HConfig header: Header event_actions: eh.EventDataAction + _last_mtime = 0. + def __init__(self, config: E2HConfig): self.config = config @@ -28,23 +38,144 @@ class E2HReduction: """ Does not do any actual reduction. """ + self.path_resolver = PathResolver(self.config.reader.year, self.config.reader.rawPath) + self.file_list = self.path_resolver.resolve(self.config.reduction.fileIdentifier) + self.file_index = 0 + self.plot_kwds = {} + self.fig = plt.figure() + + if self.config.reduction.update: + # live update implies plotting + self.config.reduction.show_plot = True + # Actions on datasets not used for normalization self.event_actions = eh.ApplyPhaseOffset(self.config.experiment.chopperPhaseOffset) - self.event_actions |= eh.CorrectChopperPhase() - self.event_actions |= ea.ExtractWalltime() + if not self.config.reduction.fast: + self.event_actions |= eh.CorrectChopperPhase() + self.event_actions |= ea.ExtractWalltime() + else: + logging.info(' Fast reduction always uses time normalization') + self.config.experiment.monitorType = MonitorType.time self.event_actions |= eh.AssociatePulseWithMonitor(self.config.experiment.monitorType) - if self.config.experiment.monitorType in [MonitorType.proton_charge or MonitorType.debug]: + if self.config.experiment.monitorType in [MonitorType.proton_charge, MonitorType.debug]: # the filtering only makes sense if using actual monitor data, not time self.event_actions |= eh.FilterMonitorThreshold(self.config.experiment.lowCurrentThreshold) - self.event_actions |= eh.FilterStrangeTimes() - self.event_actions |= ea.MergeFrames() - self.event_actions |= ea.AnalyzePixelIDs(self.config.experiment.yRange) - self.event_actions |= eh.TofTimeCorrection(self.config.experiment.incidentAngle==IncidentAngle.alphaF) - self.event_actions |= ea.CalculateWavelength(self.config.experiment.lambdaRange) + if not self.config.reduction.fast: + self.event_actions |= eh.FilterStrangeTimes() + # select needed actions in depenence of plots + if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q, + E2HPlotSelection.L]: + self.event_actions |= ea.MergeFrames() + self.event_actions |= ea.AnalyzePixelIDs(self.config.experiment.yRange) + self.event_actions |= eh.TofTimeCorrection(self.config.experiment.incidentAngle==IncidentAngle.alphaF) + self.event_actions |= ea.CalculateWavelength(self.config.experiment.lambdaRange) self.event_actions |= eh.ApplyMask() - self.grid = LZGrid([0.01], [0.0, 0.25]) + # plot dependant options + if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]: + self.grid = LZGrid(0.01, [0.0, 0.25]) + + if self.config.reduction.plot==E2HPlotSelection.LT: + self.plot_kwds['colorbar'] = True + self.plot_kwds['cmap'] = str(self.config.reduction.plot_colormap) def reduce(self): self.norm = LZNormalisation.unity(self.grid) + self.prepare_graphs() + + while self.file_index < len(self.file_list): + self.read_data() + self.add_data() + + if self.config.reduction.plotArgs==E2HPlotArguments.OutputFile: + self.create_file_output() + if self.config.reduction.plotArgs!=E2HPlotArguments.OutputFile or self.config.reduction.show_plot: + self.create_graph() + + if self.config.reduction.plotArgs==E2HPlotArguments.Default: + plt.savefig(f'e2h_{self.config.reduction.plot}.png', dpi=300) + if self.config.reduction.update: + self.timer = self.fig.canvas.new_timer(1000) + self.timer.add_callback(self.update) + self.timer.start() + if self.config.reduction.show_plot: + plt.show() + + def prepare_graphs(self): + last_file_header = AmorHeader(self.file_list[-1]) + tthh = last_file_header.geometry.nu - last_file_header.geometry.mu + + + if self.config.reduction.plot==E2HPlotSelection.LT: + self.projection = LZProjection(tthh, self.grid) + if not self.config.reduction.fast: + self.projection.correct_gravity(last_file_header.geometry.detectorDistance) + + def read_data(self): + fileName = self.file_list[self.file_index] + self.file_index += 1 + self.dataset = AmorEventData(fileName) + self.event_actions(self.dataset) + self.dataset.update_header(self.header) + + self.header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1], + timestamp=self.dataset.fileDate)) + + def add_data(self): + self.monitor = self.dataset.data.pulses.monitor.sum() + self.projection.project(self.dataset, monitor=self.monitor) + + def create_file_output(self): + ... + + def create_title(self): + output = "Events to Histogram - " + output += ",".join(["#"+os.path.basename(fi)[9:15].lstrip('0') for fi in self.file_list]) + output += f"\n$\\mu$={self.dataset.geometry.mu:.2f}" + output += f" $\\nu$={self.dataset.geometry.nu:.2f}" + if self.config.reduction.update: + output += f"\n at "+datetime.now().strftime("%m/%d/%Y %I:%M:%S") + return output + + def create_graph(self): + plt.title(self.create_title()) + self.projection.plot(**self.plot_kwds) + plt.tight_layout() + + + def replace_dataset(self, latest): + self.file_list = self.path_resolver.resolve(f'{latest}') + self.file_index = 0 + self.read_data() + self.projection.clear() + self.add_data() + self.fig.clear() + self.create_graph() + + def update(self): + if self.config.reduction.fileIdentifier=='0': + # if latest file was choosen, check if new one available and switch to it + current = int(os.path.basename(self.file_list[-1])[9:15]) + latest = self.path_resolver.search_latest(0) + if latest>current: + self.replace_dataset(latest) + return + if os.path.getmtime(self.file_list[-1])<=self._last_mtime: + return + + self._last_mtime = os.path.getmtime(self.file_list[-1]) + try: + update_data = AmorEventData(self.file_list[-1], self.dataset.last_index+1) + except EOFError: + return + + + self.event_actions(update_data) + self.dataset=update_data + self.monitor = self.dataset.data.pulses.monitor.sum() + self.projection.project(update_data, self.monitor) + + self.projection.update_plot() + plt.title(self.create_title()) + plt.draw() \ No newline at end of file diff --git a/eos/reduction_reflectivity.py b/eos/reduction_reflectivity.py index 94d316a..52860a7 100644 --- a/eos/reduction_reflectivity.py +++ b/eos/reduction_reflectivity.py @@ -158,7 +158,7 @@ class ReflectivityReduction: self.analyze_unsliced(i) def analyze_unsliced(self, i): - self.monitor = np.sum(self.dataset.data.pulses.monitor) + self.monitor = self.dataset.data.pulses.monitor.sum() logging.warning(f' monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.config.experiment.monitorType]}') proj:LZProjection = self.project_on_lz() diff --git a/setup.cfg b/setup.cfg index b69a415..7a88dc9 100644 --- a/setup.cfg +++ b/setup.cfg @@ -34,3 +34,4 @@ Homepage = "https://github.com/jochenstahn/amor" [options.entry_points] console_scripts = eos = eos.__main__:main + events2histogram = eos.e2h:main