diff --git a/libeos/event_analysis.py b/libeos/event_analysis.py index dae5894..b772991 100644 --- a/libeos/event_analysis.py +++ b/libeos/event_analysis.py @@ -14,7 +14,7 @@ from .options import IncidentAngle from .header import Header -class AnalyzedEventData(EventDataAction): +class AnalyzePixelIDs(EventDataAction): def __init__(self, yRange: Tuple[int, int]): self.yRange = yRange diff --git a/libeos/event_data_types.py b/libeos/event_data_types.py index bf58a37..bb3c381 100644 --- a/libeos/event_data_types.py +++ b/libeos/event_data_types.py @@ -1,7 +1,7 @@ """ Specify the data type and protocol used for event datasets. """ -from typing import Optional, Protocol +from typing import List, Optional, Protocol from dataclasses import dataclass from .header import Header from abc import ABC, abstractmethod @@ -60,10 +60,19 @@ class EventDatasetProtocol(Protocol): timing: AmorTiming data: AmorEventStream + def append(self, other): + # Should define a way to add events from other to own + ... + + def update_header(self, header:Header): + # update a header with the information read from file + ... + class EventDataAction(ABC): """ Abstract base class used for actions applied to an EventDatasetProtocol based objects. Each action can optionally modify the header information. + Actions can be combined using the pipe operator | (OR). """ def __call__(self, dataset: EventDatasetProtocol)->None: @@ -77,3 +86,36 @@ class EventDataAction(ABC): if hasattr(self, 'action_name'): header.reduction.corrections.append(getattr(self, 'action_name')) + def __or__(self, other:'EventDataAction')->'CombinedAction': + return CombinedAction([self, other]) + + def __repr__(self): + output = self.__class__.__name__+'(' + for key,value in self.__dict__.items(): + output += f'{key}={value}, ' + return output.rstrip(', ')+')' + +class CombinedAction(EventDataAction): + """ + Used to perform multiple actions in one call. Stores a sequence of actions + that are then performed individually one after the other. + """ + def __init__(self, actions: List[EventDataAction]) -> None: + self._actions = actions + + def perform_action(self, dataset: EventDatasetProtocol)->None: + for action in self._actions: + action(dataset) + + def update_header(self, header:Header)->None: + for action in self._actions: + action.update_header(header) + + def __or__(self, other:'EventDataAction')->'CombinedAction': + return CombinedAction(self._actions+[other]) + + def __repr__(self): + output = repr(self._actions[0]) + for ai in self._actions[1:]: + output += ' | '+repr(ai) + return output diff --git a/libeos/event_handling.py b/libeos/event_handling.py index ef4c257..6d36eda 100644 --- a/libeos/event_handling.py +++ b/libeos/event_handling.py @@ -4,6 +4,7 @@ Calculations performed on AmorEventData. import logging import numpy as np +from . import const from .options import MonitorType from .event_data_types import EventDatasetProtocol, EventDataAction from .helpers import merge_frames @@ -47,17 +48,21 @@ class AssociatePulseWithMonitor(EventDataAction): np.savetxt('tme.hst', np.vstack((dataset.data.pulses.time[:-1], cpp, dataset.data.pulses.monitor[:-1])).T) if self.monitorType in [MonitorType.proton_charge or MonitorType.debug]: - goodTimeS = dataset.data.pulses.time[dataset.data.pulses.monitor!=0] - filter_e = np.isin(dataset.data.events.wallTime, goodTimeS) - dataset.data.events = dataset.data.events[filter_e] - logging.info(f' low-beam (<{self.lowCurrentThreshold} mC) rejected pulses: ' - f'{dataset.data.pulses.monitor.shape[0]-goodTimeS.shape[0]} ' - f'out of {dataset.data.pulses.monitor.shape[0]}') - logging.info(f' with {filter_e.shape[0]-dataset.data.events.shape[0]} events') - if goodTimeS.shape[0]: - logging.info(f' average counts per pulse = {dataset.data.events.shape[0] / goodTimeS.shape[0]:7.1f}') - else: - logging.info(f' average counts per pulse = undefined') + self.monitor_threshold(dataset) + + def monitor_threshold(self, dataset): + # TODO: evaluate if this should actually do masking instead + goodTimeS = dataset.data.pulses.time[dataset.data.pulses.monitor!=0] + filter_e = np.isin(dataset.data.events.wallTime, goodTimeS) + dataset.data.events = dataset.data.events[filter_e] + logging.info(f' low-beam (<{self.lowCurrentThreshold} mC) rejected pulses: ' + f'{dataset.data.pulses.monitor.shape[0]-goodTimeS.shape[0]} ' + f'out of {dataset.data.pulses.monitor.shape[0]}') + logging.info(f' with {filter_e.shape[0]-dataset.data.events.shape[0]} events') + if goodTimeS.shape[0]: + logging.info(f' average counts per pulse = {dataset.data.events.shape[0]/goodTimeS.shape[0]:7.1f}') + else: + logging.info(f' average counts per pulse = undefined') @staticmethod def get_current_per_pulse(pulseTimeS, currentTimeS, currents): @@ -81,10 +86,8 @@ class FilterStrangeTimes(EventDataAction): logging.warning(f' strange times: {np.logical_not(filter_e).sum()}') class MergeFrames(EventDataAction): - def __init__(self, tofCut:float): - self.tofCut = tofCut - def perform_action(self, dataset: EventDatasetProtocol)->None: - total_offset = (self.tofCut + + tofCut = const.lamdaCut+dataset.geometry.chopperDetectorDistance/const.hdm*1e-13 + total_offset = (tofCut + dataset.timing.tau * (dataset.timing.ch1TriggerPhase + dataset.timing.chopperPhase/2)/180) dataset.data.events.tof = merge_frames(dataset.data.events.tof, self.tofCut, dataset.timing.tau, total_offset) diff --git a/libeos/file_reader.py b/libeos/file_reader.py index b024d73..7d7c995 100644 --- a/libeos/file_reader.py +++ b/libeos/file_reader.py @@ -1,23 +1,26 @@ """ Reading of Amor NeXus data files to extract metadata and event stream. """ -from typing import BinaryIO, Union +from typing import BinaryIO, List, Union import h5py import numpy as np import platform import logging import subprocess +import sys +import os from datetime import datetime from orsopy import fileio from orsopy.fileio.model_language import SampleModel -from . import const +from . import const, event_handling as eh, event_analysis as ea from .header import Header from .helpers import extract_walltime from .event_data_types import AmorGeometry, AmorTiming, AmorEventStream, PACKET_TYPE, EVENT_TYPE, PULSE_TYPE, PC_TYPE +from .options import ExperimentConfig, IncidentAngle, ReaderConfig try: import zoneinfo @@ -35,6 +38,47 @@ if platform.node().startswith('amor'): else: NICOS_CACHE_DIR = None +class PathResolver: + def __init__(self, year, rawPath): + self.year = year + self.rawPath = rawPath + + def resolve(self, short_notation): + return list(map(self.get_path, self.expand_file_list(short_notation))) + + def expand_file_list(self, short_notation): + """Evaluate string entry for file number lists""" + file_list = [] + for i in short_notation.split(','): + if '-' in i: + if ':' in i: + step = i.split(':', 1)[1] + file_list += range(int(i.split('-', 1)[0]), + int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, + int(step)) + else: + step = 1 + file_list += range(int(i.split('-', 1)[0]), + int(i.split('-', 1)[1])+1, + int(step)) + else: + file_list += [int(i)] + return list(sorted(file_list)) + + def get_path(self, number): + fileName = f'amor{self.year}n{number:06d}.hdf' + path = '' + for rawd in self.rawPath: + if os.path.exists(os.path.join(rawd, fileName)): + 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)}' + else: + sys.exit(f'# ERROR: the file {fileName} can not be found in {self.rawPath}') + return os.path.join(path, fileName) class AmorEventData: """ @@ -55,7 +99,7 @@ class AmorEventData: timing: AmorTiming data: AmorEventStream - startTime: np.int64 + eventStartTime: np.int64 def __init__(self, fileName:Union[str, h5py.File, BinaryIO], first_index:int=0, max_events:int=1_000_000): if type(fileName) is str: @@ -301,8 +345,119 @@ class AmorEventData: output += '\n' return output + def append(self, other): + """ + Append event streams from another file to this one. Adjusts the event indices in the + packets to stay valid. + """ + new_events = np.concatenate([self.data.events, other.data.events]).view(np.recarray) + new_pulses = np.concatenate([self.data.pulses, other.data.pulses]).view(np.recarray) + new_proton_current = np.concatenate([self.data.proton_current, other.data.proton_current]).view(np.recarray) + new_packets = np.concatenate([self.data.packets, other.data.packets]).view(np.recarray) + new_packets.start_index[self.data.packets.shape[0]:] += self.data.events.shape[0] + self.data = AmorEventStream(new_events, new_packets, new_pulses, new_proton_current) + # Indicate that this is amodified dataset, basically counts number of appends as negative indices + self.last_index = min(self.last_index-1, -1) + + def __repr__(self): output = (f"AmorEventData({self.fileName!r}) # {self.data.events.shape[0]} events, " f"{self.data.pulses.shape[0]} pulses") return output + +class AmorData: + """re-implement old AmorData class functionality until refactoring is complete""" + chopperDetectorDistance: float + chopperDistance: float + chopperPhase: float + chopperSpeed: float + chopper1TriggerPhase: float + chopper2TriggerPhase: float + div: float + data_file_numbers: List[int] + delta_z: np.ndarray + detZ_e: np.ndarray + lamda_e: np.ndarray + wallTime_e: np.ndarray + kad: float + kap: float + lambdaMax: float + lambda_e: np.ndarray + # monitor: float + mu: float + nu: float + tau: float + tofCut: float + start_date: str + monitorType: str + + seriesStartTime = None + + # ------------------------------------------------------------------------------------------------- + def __init__(self, header: Header, reader_config: ReaderConfig, config: ExperimentConfig, + short_notation: str, norm=False): + # self.startTime = reader_config.startTime + self.header = header + self.config = config + self.reader_config = reader_config + resolver = PathResolver(reader_config.year, reader_config.rawPath) + self.file_list = resolver.resolve(short_notation) + self.norm = norm + self.prepare_actions() + self.process() + self.assign() + + def prepare_actions(self): + # setup all actions performed in origin AmorData, time correction requires first dataset start time + self.event_actions = eh.AssociatePulseWithMonitor(self.config.monitorType, self.config.lowCurrentThreshold) + self.event_actions |= eh.FilterStrangeTimes() + self.event_actions |= eh.MergeFrames() + self.event_actions |= ea.AnalyzePixelIDs(self.config.yRange) + self.event_actions |= ea.TofTimeCorrection(self.config.incidentAngle==IncidentAngle.alphaF) + self.event_actions |= ea.WavelengthAndQ(self.config.lambdaRange, self.config.incidentAngle) + self.event_actions |= ea.FilterQzRange(self.config.qzRange) + self.event_actions |= ea.ApplyMask() + + def process(self): + self.dataset = AmorEventData(self.file_list[0]) + time_correction = eh.CorrectSeriesTime(self.dataset.eventStartTime) + time_correction(self.dataset) + self.event_actions(self.dataset) + if not self.norm: + self.dataset.update_header(self.header) + time_correction.update_header(self.header) + self.event_actions.update_header(self.header) + for fi in self.file_list[1:]: + di = AmorEventData(fi) + time_correction(di) + self.event_actions(di) + self.dataset.append(di) + + for fileName in self.file_list: + if self.norm: + self.header.measurement_additional_files.append(fileio.File( + file=fileName.split('/')[-1], + timestamp=self.dataset.fileDate)) + else: + self.header.measurement_data_files.append(fileio.File( + file=fileName.split('/')[-1], + timestamp=self.dataset.fileDate)) + + def assign(self): + # assigne old class parameters from dataset object. + ds = self.dataset + ev = ds.data.events + self.detZ_e = ev.detZ + self.lamda_e = ev.lamda + self.wallTime_e = ev.wallTime + self.qz_e = ev.qz + self.qx_e = ev.qx + self.pulseTimeS = ds.data.pulses.time + self.monitorPerPulse = ds.data.pulses.monitor + + for key, value in ds.geometry.__dict__.items(): + setattr(self, key, value) + for key, value in ds.timing.__dict__.items(): + setattr(self, key, value) + self.startTime = ds.eventStartTime \ No newline at end of file diff --git a/libeos/options.py b/libeos/options.py index e98bead..4164d8b 100644 --- a/libeos/options.py +++ b/libeos/options.py @@ -197,8 +197,8 @@ class ExperimentConfig(ArgParsable): 'help': 'rotation speed of the chopper disks in rpm', }, ) - yRange: Tuple[float, float] = field( - default_factory=lambda: [18, 48], + yRange: Tuple[int, int] = field( + default=(18, 48), metadata={ 'short': 'y', 'group': 'region of interest', diff --git a/libeos/reduction.py b/libeos/reduction.py index cf37ed0..e0987b4 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -5,7 +5,7 @@ import sys import numpy as np from orsopy import fileio -from .file_reader_old import AmorData +from .file_reader import AmorData from .header import Header from .options import EOSConfig, IncidentAngle, MonitorType, NormalisationMethod from .instrument import Grid