diff --git a/libeos/file_reader_old.py b/libeos/file_reader_old.py index 820f60f..f7ec7c0 100644 --- a/libeos/file_reader_old.py +++ b/libeos/file_reader_old.py @@ -2,17 +2,17 @@ import logging import os import subprocess import sys -import platform from datetime import datetime, timezone -from dataclasses import dataclass + + try: import zoneinfo except ImportError: # for python versions < 3.9 try to use the backports version from backports import zoneinfo -from typing import List, Optional -from abc import ABC, abstractmethod +from typing import List +import yaml import h5py import numpy as np from orsopy import fileio @@ -22,16 +22,16 @@ from . import const from .header import Header from .instrument import Detector from .options import ExperimentConfig, IncidentAngle, MonitorType, ReaderConfig -from .helpers import merge_frames, extract_walltime, filter_project_x, calculate_derived_properties_focussing + + +try: + from . import helpers_numba as nb_helpers +except Exception: + nb_helpers = None # Time zone used to interpret time strings 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"' -else: - NICOS_CACHE_DIR = None class AmorData: """read meta-data and event streams from .hdf file(s), apply filters and conversions""" @@ -51,7 +51,7 @@ class AmorData: kap: float lambdaMax: float lambda_e: np.ndarray - #monitor: float + # monitor: float mu: float nu: float tau: float @@ -61,17 +61,17 @@ class AmorData: seriesStartTime = None - #------------------------------------------------------------------------------------------------- + # ------------------------------------------------------------------------------------------------- def __init__(self, header: Header, reader_config: ReaderConfig, config: ExperimentConfig, - short_notation:str, norm=False): - #self.startTime = reader_config.startTime + short_notation: str, norm=False): + # self.startTime = reader_config.startTime self.header = header self.config = config self.reader_config = reader_config self.expand_file_list(short_notation) self.read_data(norm=norm) - #------------------------------------------------------------------------------------------------- + # ------------------------------------------------------------------------------------------------- def read_data(self, norm=False): self.file_list = [] for number in self.data_file_numbers: @@ -85,7 +85,7 @@ class AmorData: _detZ_e = [] _lamda_e = [] _wallTime_e = [] - #_monitor = 0 + # _monitor = 0 _monitorPerPulse = [] _pulseTimeS = [] for file in self.file_list: @@ -95,33 +95,35 @@ class AmorData: _wallTime_e = np.append(_wallTime_e, self.wallTime_e) _monitorPerPulse = np.append(_monitorPerPulse, self.monitorPerPulse) _pulseTimeS = np.append(_pulseTimeS, self.pulseTimeS) - #_monitor += self.monitor + # _monitor += self.monitor self.detZ_e = _detZ_e self.lamda_e = _lamda_e self.wallTime_e = _wallTime_e - #self.monitor = _monitor - self.monitorPerPulse = _monitorPerPulse - self.pulseTimeS = _pulseTimeS + # self.monitor = _monitor + self.monitorPerPulse = _monitorPerPulse + self.pulseTimeS = _pulseTimeS - #------------------------------------------------------------------------------------------------- + # ------------------------------------------------------------------------------------------------- def path_generator(self, number): fileName = f'amor{self.reader_config.year}n{number:06d}.hdf' path = '' for rawd in self.reader_config.rawPath: - if os.path.exists(os.path.join(rawd,fileName)): + 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.reader_config.year}/amor/{int(number/1000)}/{fileName}'): + if os.path.exists( + f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}/{fileName}'): path = f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}' else: sys.exit(f'# ERROR: the file {fileName} can not be found in {self.reader_config.rawPath}') return os.path.join(path, fileName) - #------------------------------------------------------------------------------------------------- + + # ------------------------------------------------------------------------------------------------- def expand_file_list(self, short_notation): """Evaluate string entry for file number lists""" - #log().debug('Executing get_flist') - file_list=[] + # log().debug('Executing get_flist') + file_list = [] for i in short_notation.split(','): if '-' in i: if ':' in i: @@ -136,23 +138,25 @@ class AmorData: int(step)) else: file_list += [int(i)] - self.data_file_numbers=sorted(file_list) - #------------------------------------------------------------------------------------------------- + self.data_file_numbers = sorted(file_list) + + # ------------------------------------------------------------------------------------------------- def resolve_pixels(self): """determine spatial coordinats and angles from pixel number""" - nPixel = Detector.nWires * Detector.nStripes * Detector.nBlades + nPixel = Detector.nWires*Detector.nStripes*Detector.nBlades pixelID = np.arange(nPixel) - (bladeNr, bPixel) = np.divmod(pixelID, Detector.nWires * Detector.nStripes) - (bZi, detYi) = np.divmod(bPixel, Detector.nStripes) # z index on blade, y index on detector - detZi = bladeNr * Detector.nWires + bZi # z index on detector - detX = bZi * Detector.dX # x position in detector + (bladeNr, bPixel) = np.divmod(pixelID, Detector.nWires*Detector.nStripes) + (bZi, detYi) = np.divmod(bPixel, Detector.nStripes) # z index on blade, y index on detector + detZi = bladeNr*Detector.nWires+bZi # z index on detector + detX = bZi*Detector.dX # x position in detector # detZ = Detector.zero - bladeNr * Detector.bladeZ - bZi * Detector.dZ # z position on detector - bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*Detector.bladeZ / Detector.distance) ) - delta = (Detector.nBlades/2. - bladeNr) * bladeAngle \ - - np.rad2deg( np.arctan(bZi*Detector.dZ / ( Detector.distance + bZi * Detector.dX) ) ) - self.delta_z = delta[detYi==1] + bladeAngle = np.rad2deg(2.*np.arcsin(0.5*Detector.bladeZ/Detector.distance)) + delta = (Detector.nBlades/2.-bladeNr)*bladeAngle \ + -np.rad2deg(np.arctan(bZi*Detector.dZ/(Detector.distance+bZi*Detector.dX))) + self.delta_z = delta[detYi==1] return np.vstack((detYi.T, detZi.T, detX.T, delta.T)).T - #------------------------------------------------------------------------------------------------- + + # ------------------------------------------------------------------------------------------------- def read_individual_data(self, fileName, norm=False): self.hdf = h5py.File(fileName, 'r', swmr=True) @@ -166,13 +170,13 @@ class AmorData: if self.readHeaderInfo: self.readHeaderInfo = False self.header.measurement_instrument_settings = fileio.InstrumentSettings( - incident_angle = fileio.ValueRange(round(self.mu+self.kap+self.kad-0.5*self.div, 3), - round(self.mu+self.kap+self.kad+0.5*self.div, 3), - 'deg'), - wavelength = fileio.ValueRange(const.lamdaCut, self.config.lambdaRange[1], 'angstrom'), - #polarization = fileio.Polarization.unpolarized, - polarization = fileio.Polarization(self.polarizationConfig) - ) + incident_angle=fileio.ValueRange(round(self.mu+self.kap+self.kad-0.5*self.div, 3), + round(self.mu+self.kap+self.kad+0.5*self.div, 3), + 'deg'), + wavelength=fileio.ValueRange(const.lamdaCut, self.config.lambdaRange[1], 'angstrom'), + # polarization = fileio.Polarization.unpolarized, + polarization=fileio.Polarization(self.polarizationConfig) + ) self.header.measurement_instrument_settings.mu = fileio.Value( round(self.mu, 3), 'deg', @@ -196,12 +200,12 @@ class AmorData: comment='incoming beam angular offset') if norm: self.header.measurement_additional_files.append(fileio.File( - file=fileName.split('/')[-1], - timestamp=self.fileDate)) + file=fileName.split('/')[-1], + timestamp=self.fileDate)) else: self.header.measurement_data_files.append(fileio.File( - file=fileName.split('/')[-1], - timestamp=self.fileDate)) + file=fileName.split('/')[-1], + timestamp=self.fileDate)) logging.info(f' mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kad:6.3f}') @@ -218,11 +222,11 @@ class AmorData: self.associate_pulse_with_monitor() # following lines: debugging output to trace the time-offset of proton current and neutron pulses - if self.config.monitorType == MonitorType.debug: + if self.config.monitorType==MonitorType.debug: cpp, t_bins = np.histogram(self.wallTime_e, self.pulseTimeS) np.savetxt('tme.hst', np.vstack((self.pulseTimeS[:-1], cpp, self.monitorPerPulse[:-1])).T) - #self.average_events_per_pulse() # for debugging only. VERY time consuming!!! + # self.average_events_per_pulse() # for debugging only. VERY time consuming!!! self.monitor_threshold() @@ -238,7 +242,8 @@ class AmorData: self.filter_qz_range(norm) - logging.info(f' number of events: total = {self.totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') + logging.info( + f' number of events: total = {self.totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') def read_event_stream(self): self.tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 @@ -247,15 +252,15 @@ class AmorData: self.dataPacketTime_p = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.int64) def correct_for_chopper_phases(self): - #print(f'tof phase-offset: {self.ch1TriggerPhase - self.chopperPhase/2}') - self.tof_e += self.tau * (self.ch1TriggerPhase - self.chopperPhase/2)/180 + # print(f'tof phase-offset: {self.ch1TriggerPhase - self.chopperPhase/2}') + self.tof_e += self.tau*(self.ch1TriggerPhase-self.chopperPhase/2)/180 def read_chopper_trigger_stream(self): self.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) + # 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) - if np.shape(self.chopper1TriggerTime)[0] > 2: + if np.shape(self.chopper1TriggerTime)[0]>2: self.startTime = self.chopper1TriggerTime[0] self.stopTime = self.chopper1TriggerTime[-1] self.pulseTimeS = self.chopper1TriggerTime @@ -263,7 +268,7 @@ class AmorData: logging.warn(' no chopper trigger data available, using event steram instead') self.startTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][0], dtype=np.int64) self.stopTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-2], dtype=np.int64) - self.pulseTimeS = np.arange(self.startTime, self.stopTime, self.tau*1e9) + self.pulseTimeS = np.arange(self.startTime, self.stopTime, self.tau*1e9) if self.seriesStartTime is None: self.seriesStartTime = self.startTime logging.debug(f' series start time (epoch): {self.seriesStartTime/1e9:13.2f} s') @@ -272,15 +277,21 @@ class AmorData: logging.debug(f' => counting time {self.stopTime/1e9-self.startTime/1e9:8.2f} s') def extract_walltime(self, norm): - self.wallTime_e = extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p) + if nb_helpers: + self.wallTime_e = nb_helpers.extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p) + else: + self.wallTime_e = np.empty(np.shape(self.tof_e)[0], dtype=np.int64) + for i in range(len(self.dataPacket_p)-1): + self.wallTime_e[self.dataPacket_p[i]:self.dataPacket_p[i+1]] = self.dataPacketTime_p[i] + self.wallTime_e[self.dataPacket_p[-1]:] = self.dataPacketTime_p[-1] self.wallTime_e -= np.int64(self.seriesStartTime) logging.debug(f' wall time from {self.wallTime_e[0]/1e9:6.1f} s to {self.wallTime_e[-1]/1e9:6.1f} s') def read_proton_current_stream(self): self.currentTime = np.array(self.hdf['entry1/Amor/detector/proton_current/time'][:], dtype=np.int64) - self.current = np.array(self.hdf['entry1/Amor/detector/proton_current/value'][:,0], dtype=float) - if self.config.monitorType == MonitorType.auto: - if self.current.sum() > 1: + self.current = np.array(self.hdf['entry1/Amor/detector/proton_current/value'][:, 0], dtype=float) + if self.config.monitorType==MonitorType.auto: + if self.current.sum()>1: self.monitorType = MonitorType.proton_charge logging.warn(' monitor type set to "proton current"') else: @@ -288,19 +299,19 @@ class AmorData: logging.warn(' monitor type set to "time"') def associate_pulse_with_monitor(self): - if self.config.monitorType == MonitorType.proton_charge or MonitorType.debug: + if self.config.monitorType==MonitorType.proton_charge or MonitorType.debug: self.currentTime -= np.int64(self.seriesStartTime) self.monitorPerPulse = self.get_current_per_pulse(self.pulseTimeS, self.currentTime, - self.current)\ - * 2*self.tau * 1e-3 + self.current) \ + *2*self.tau*1e-3 # filter low-current pulses - self.monitorPerPulse = np.where(self.monitorPerPulse > 2*self.tau * self.config.lowCurrentThreshold * 1e-3, + self.monitorPerPulse = np.where(self.monitorPerPulse>2*self.tau*self.config.lowCurrentThreshold*1e-3, self.monitorPerPulse, 0) - elif self.config.monitorType == MonitorType.time: + elif self.config.monitorType==MonitorType.time: self.monitorPerPulse = np.ones(np.shape(self.pulseTimeS)[0])*2*self.tau - else: # pulses + else: # pulses self.monitorPerPulse = np.ones(np.shape(self.pulseTimeS)[0]) def get_current_per_pulse(self, pulseTimeS, currentTimeS, currents): @@ -310,19 +321,19 @@ class AmorData: pulseCurrentS = np.zeros(pulseTimeS.shape[0], dtype=float) j = 0 for i, ti in enumerate(pulseTimeS): - while ti >= currentTimeS[j+1]: + while ti>=currentTimeS[j+1]: j += 1 pulseCurrentS[i] = currents[j] return pulseCurrentS def average_events_per_pulse(self): - if self.config.monitorType == MonitorType.proton_charge: + if self.config.monitorType==MonitorType.proton_charge: for i, time in enumerate(self.pulseTimeS): - events = np.shape(self.wallTime_e[self.wallTime_e == time])[0] + events = np.shape(self.wallTime_e[self.wallTime_e==time])[0] logging.info(f'pulse: {i:6.0f}, events: {events:6.0f}, monitor: {self.monitorPerPulse[i]:6.2f}') def monitor_threshold(self): - #if self.config.monitorType == MonitorType.proton_charge: # fix to check for file compatibility + # if self.config.monitorType == MonitorType.proton_charge: # fix to check for file compatibility self.totalNumber = np.shape(self.tof_e[self.tof_e<=self.stopTime])[0] if True: goodTimeS = self.pulseTimeS[self.monitorPerPulse!=0] @@ -330,10 +341,12 @@ class AmorData: self.tof_e = self.tof_e[filter_e] self.pixelID_e = self.pixelID_e[filter_e] self.wallTime_e = self.wallTime_e[filter_e] - logging.info(f' low-beam (<{self.config.lowCurrentThreshold} mC) rejected pulses: {np.shape(self.monitorPerPulse)[0]-1-np.shape(goodTimeS)[0]} out of {np.shape(self.monitorPerPulse)[0]-1}') + logging.info( + f' low-beam (<{self.config.lowCurrentThreshold} mC) rejected pulses: {np.shape(self.monitorPerPulse)[0]-1-np.shape(goodTimeS)[0]} out of {np.shape(self.monitorPerPulse)[0]-1}') logging.info(f' with {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]} events') if np.shape(goodTimeS[goodTimeS!=0])[0]: - logging.info(f' average counts per pulse = {np.shape(self.tof_e)[0] / np.shape(goodTimeS[goodTimeS!=0])[0]:7.1f}') + logging.info( + f' average counts per pulse = {np.shape(self.tof_e)[0]/np.shape(goodTimeS[goodTimeS!=0])[0]:7.1f}') else: logging.info(f' average counts per pulse = undefined') @@ -348,27 +361,38 @@ class AmorData: # TODO: - handle each neutron pulse individually, - associate with correct monitor also for slow neutrons def merge_time_frames(self): - total_offset = self.tofCut + self.tau * (self.ch1TriggerPhase + self.chopperPhase/2)/180 - self.tof_e = merge_frames(self.tof_e, self.tofCut, self.tau, total_offset) + total_offset = self.tofCut+self.tau*(self.ch1TriggerPhase+self.chopperPhase/2)/180 + if nb_helpers: + self.tof_e = nb_helpers.merge_frames(self.tof_e, self.tofCut, self.tau, total_offset) + else: + self.tof_e = np.remainder(self.tof_e-(self.tofCut-self.tau), + self.tau)+total_offset # tof shifted to 1 frame def filter_project_x(self): pixelLookUp = self.resolve_pixels() - (self.detZ_e, self.detXdist_e, self.delta_e, self.mask_e) = filter_project_x( - pixelLookUp, self.pixelID_e.astype(np.int64), self.config.yRange[0], self.config.yRange[1] - ) + if nb_helpers: + (self.detZ_e, self.detXdist_e, self.delta_e, self.mask_e) = nb_helpers.filter_project_x( + pixelLookUp, self.pixelID_e.astype(np.int64), self.config.yRange[0], self.config.yRange[1] + ) + else: + # resolve pixel ID into y and z indicees, x position and angle + (detY_e, self.detZ_e, self.detXdist_e, self.delta_e) = pixelLookUp[np.int_(self.pixelID_e)-1, :].T + # define mask and filter y range + self.mask_e = (self.config.yRange[0]<=detY_e) & (detY_e<=self.config.yRange[1]) def correct_for_chopper_opening(self): # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau - if self.config.incidentAngle == IncidentAngle.alphaF: - self.tof_e -= ( self.delta_e / 180. ) * self.tau + if self.config.incidentAngle==IncidentAngle.alphaF: + self.tof_e -= (self.delta_e/180.)*self.tau else: # TODO: check sign of correction - self.tof_e -= ( self.kad / 180. ) * self.tau + self.tof_e -= (self.kad/180.)*self.tau def calculate_derived_properties(self): self.lamdaMax = const.lamdaCut+1.e13*self.tau*const.hdm/(self.chopperDetectorDistance+124.) + # if nb_helpers: if False: - self.lamda_e, self.qz_e, self.mask_e = calculate_derived_properties_focussing( + self.lamda_e, self.qz_e, self.mask_e = nb_helpers.calculate_derived_properties_focussing( self.tof_e, self.detXdist_e, self.delta_e, self.mask_e, self.config.lambdaRange[0], self.config.lambdaRange[1], self.nu, self.mu, self.chopperDetectorDistance, const.hdm @@ -377,24 +401,24 @@ class AmorData: # lambda self.lamda_e = (1.e13*const.hdm)*self.tof_e/(self.chopperDetectorDistance+self.detXdist_e) self.mask_e = np.logical_and(self.mask_e, (self.config.lambdaRange[0]<=self.lamda_e) & ( - self.lamda_e<=self.config.lambdaRange[1])) + self.lamda_e<=self.config.lambdaRange[1])) # alpha_f # q_z - if self.config.incidentAngle == IncidentAngle.alphaF: - alphaF_e = self.nu - self.mu + self.delta_e + if self.config.incidentAngle==IncidentAngle.alphaF: + alphaF_e = self.nu-self.mu+self.delta_e self.qz_e = 4*np.pi*(np.sin(np.deg2rad(alphaF_e))/self.lamda_e) # qx_e = 0. self.header.measurement_scheme = 'angle- and energy-dispersive' - elif self.config.incidentAngle == IncidentAngle.nu: - alphaF_e = (self.nu + self.delta_e + self.kap + self.kad) / 2. + elif self.config.incidentAngle==IncidentAngle.nu: + alphaF_e = (self.nu+self.delta_e+self.kap+self.kad)/2. self.qz_e = 4*np.pi*(np.sin(np.deg2rad(alphaF_e))/self.lamda_e) # qx_e = 0. self.header.measurement_scheme = 'energy-dispersive' else: - alphaF_e = self.nu - self.mu + self.delta_e - alphaI = self.kap + self.kad + self.mu - self.qz_e = 2*np.pi * ((np.sin(np.deg2rad(alphaF_e)) + np.sin(np.deg2rad(alphaI)))/self.lamda_e) - self.qx_e = 2*np.pi * ((np.cos(np.deg2rad(alphaF_e)) - np.cos(np.deg2rad(alphaI)))/self.lamda_e) + alphaF_e = self.nu-self.mu+self.delta_e + alphaI = self.kap+self.kad+self.mu + self.qz_e = 2*np.pi*((np.sin(np.deg2rad(alphaF_e))+np.sin(np.deg2rad(alphaI)))/self.lamda_e) + self.qx_e = 2*np.pi*((np.cos(np.deg2rad(alphaF_e))-np.cos(np.deg2rad(alphaI)))/self.lamda_e) self.header.measurement_scheme = 'energy-dispersive' def filter_qz_range(self, norm): @@ -405,37 +429,35 @@ class AmorData: self.lamda_e = self.lamda_e[self.mask_e] self.wallTime_e = self.wallTime_e[self.mask_e] - - def read_individual_header(self): self.chopperDistance = float(np.take(self.hdf['entry1/Amor/chopper/pair_separation'], 0)) self.detectorDistance = float(np.take(self.hdf['entry1/Amor/detector/transformation/distance'], 0)) self.chopperDetectorDistance = self.detectorDistance-float(np.take(self.hdf['entry1/Amor/chopper/distance'], 0)) self.tofCut = const.lamdaCut*self.chopperDetectorDistance/const.hdm*1.e-13 - #TODO: 'undefined' is not orso compatible - but should be. + # TODO: 'undefined' is not orso compatible - but should be. polarizationConfigs = ['undefined', 'unpolarized', 'po', 'mo', 'op', 'pp', 'mp', 'om', 'pm', 'mm'] try: - self.mu = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/mu'], 0)) - self.nu = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/nu'], 0)) - self.kap = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kappa'], 0)) - self.kad = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kappa_offset'], 0)) - #self.kap = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kap'], 0)) - #self.kad = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kad'], 0)) - self.div = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/div'], 0)) + self.mu = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/mu'], 0)) + self.nu = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/nu'], 0)) + self.kap = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kappa'], 0)) + self.kad = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kappa_offset'], 0)) + # self.kap = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kap'], 0)) + # self.kad = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kad'], 0)) + self.div = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/div'], 0)) self.ch1TriggerPhase = float(np.take(self.hdf['/entry1/Amor/chopper/ch1_trigger_phase'], 0)) self.ch2TriggerPhase = float(np.take(self.hdf['/entry1/Amor/chopper/ch2_trigger_phase'], 0)) - try: - chopperTriggerTime = (float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][7])\ - - float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][0]))\ - / 7 + try: + chopperTriggerTime = (float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][7]) \ + -float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][0])) \ + /7 self.tau = int(1e-6*chopperTriggerTime/2+0.5)*(1e-3) self.chopperSpeed = 30/self.tau - chopperTriggerTimeDiff = float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][2]) + chopperTriggerTimeDiff = float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][2]) chopperTriggerPhase = 180e-9*chopperTriggerTimeDiff/self.tau - #TODO: check the next line - self.chopperPhase = chopperTriggerPhase + self.ch1TriggerPhase - self.ch2TriggerPhase - #print(f'chopperTriggerPhase: {chopperTriggerPhase} + {self.ch1TriggerPhase} - {self.ch2TriggerPhase} chopper phase: {self.chopperPhase}') + # TODO: check the next line + self.chopperPhase = chopperTriggerPhase+self.ch1TriggerPhase-self.ch2TriggerPhase + # print(f'chopperTriggerPhase: {chopperTriggerPhase} + {self.ch1TriggerPhase} - {self.ch2TriggerPhase} chopper phase: {self.chopperPhase}') except(KeyError, IndexError): logging.debug(' chopper speed and phase taken from .hdf file') self.chopperSpeed = float(np.take(self.hdf['/entry1/Amor/chopper/rotation_speed'], 0)) @@ -469,10 +491,11 @@ class AmorData: self.ch1TriggerPhase = float(value) value = str(subprocess.getoutput(f'{grp} {cachePath}nicos-ch2_trigger_phase/{year_date}')).split('\t')[-1] self.ch2TriggerPhase = float(value) - value = str(subprocess.getoutput(f'{grp} {cachePath}nicos-polarizer_config_label/{year_date}')).split('\t')[-1] - self. polarizationConfigLabel = int(value) - - self.tau = 30. / self.chopperSpeed + value = str(subprocess.getoutput(f'{grp} {cachePath}nicos-polarizer_config_label/{year_date}')).split('\t')[ + -1] + self.polarizationConfigLabel = int(value) + + self.tau = 30./self.chopperSpeed self.polarizationConfig = polarizationConfigs[polarizationConfigLabel] logging.debug(f' polarization configuration: {self.polarizationConfig} (index {polarizationConfigLabel})') @@ -488,15 +511,16 @@ class AmorData: logging.debug(f' replaced nu = {self.nu} with {self.config.nu}') self.nu = self.config.nu if self.config.chopperPhaseOffset: - logging.debug(f' replaced ch1TriggerPhase = {self.ch1TriggerPhase} with {self.config.chopperPhaseOffset}') + logging.debug( + f' replaced ch1TriggerPhase = {self.ch1TriggerPhase} with {self.config.chopperPhaseOffset}') self.ch1TriggerPhase = self.config.chopperPhaseOffset # extract start time as unix time, adding UTC offset of 1h to time string dz = datetime.fromisoformat(self.hdf['/entry1/start_time'][0].decode('utf-8')) - self.fileDate=dz.replace(tzinfo=AMOR_LOCAL_TIMEZONE) - #self.startTime = np.int64( (self.fileDate.timestamp() ) * 1e9 ) - #if self.seriesStartTime is None: - # self.seriesStartTime = self.startTime + self.fileDate = dz.replace(tzinfo=AMOR_LOCAL_TIMEZONE) + # self.startTime = np.int64( (self.fileDate.timestamp() ) * 1e9 ) + # if self.seriesStartTime is None: + # self.seriesStartTime = self.startTime def read_header_info(self): # read general information and first data set @@ -509,14 +533,30 @@ class AmorData: user_email = self.hdf['entry1/user/email'][0].decode('utf-8') user_orcid = None sampleName = self.hdf['entry1/sample/name'][0].decode('utf-8') - model = self.hdf['entry1/sample/model'][0].decode('utf-8') instrumentName = 'Amor' source = self.hdf['entry1/Amor/source/name'][0].decode('utf-8') sourceProbe = 'neutron' start_time = self.hdf['entry1/start_time'][0].decode('utf-8') self.start_date = start_time.split(' ')[0] + if self.config.sampleModel: - model = self.config.sampleModel + if 'yml' in self.config.sampleModel or 'yaml' in self.config.sampleModel: + if os.path.isfile(self.config.sampleModel): + with open(self.config.sampleModel, 'r') as model_yml: + model = yaml.safe_load(model_yml) + else: + logging.warning(f' ! the file {self.config.sampleModel}.yml does not exist. Ignored!') + else: + model = dict(stack=self.config.sampleModel) + try: + model + except NameError: + _model = self.hdf['entry1/sample/model'][0].decode('utf-8') + if type(_model)==dict: + model = yaml.safe_load(_model) + else: + model = dict(stack=_model) + # assembling orso header information self.header.owner = fileio.Person( name=user_name, @@ -535,8 +575,7 @@ class AmorData: ) self.header.sample = fileio.Sample( name=sampleName, - model=SampleModel(stack=model), + model=SampleModel.from_dict(model), sample_parameters=None, ) self.header.measurement_scheme = 'angle- and energy-dispersive' - diff --git a/tests/mock_data.py b/tests/mock_data.py new file mode 100644 index 0000000..8b5f00a --- /dev/null +++ b/tests/mock_data.py @@ -0,0 +1,64 @@ +""" +Generate a mock dataset in memory for running unit tests. +""" + +import h5py +import numpy as np + +MOCK_METADATA = { + 'title': 'Testdata', + 'proposal_id': 'none', + 'user/name': 'test user', + 'user/email': 'test@user.de', + 'sample/name': 'test sample', + 'sample/model': 'air | Fe 12 | Si', + 'Amor/source/name': 'SINQ', + 'start_time': '2025-01-01 00:00:01', + } +MOCK_META_TYPED = { + 'Amor/chopper/pair_separation': (1000.0, np.float32), + 'Amor/detector/transformation/distance': (4000.0, np.float64), + 'Amor/instrument_control_parameters/kappa': (1000.0, np.float64), + 'Amor/instrument_control_parameters/kappa_offset': (1000.0, np.float64), + 'Amor/instrument_control_parameters/div': (1.6, np.float64), + 'Amor/chopper/ch1_trigger_phase': (-9.1, np.float64), + 'Amor/chopper/ch2_trigger_phase': (6.75, np.float64), + 'Amor/chopper/ch2_trigger/event_time_zero': ([0.0]*10, np.uint64), + 'Amor/chopper/ch2_trigger/event_time_offset': ([0.0]*10, np.uint32), + 'Amor/chopper/rotation_speed': (500.0, np.float64), + 'Amor/chopper/phase': (0.0, np.float64), + 'Amor/polarization/configuration/value': (0.0, np.float64), + } + +def mock_data(mu=1.0, nu=2.0): + hdf = h5py.File.in_memory() # requires h5py >=3.13 + ds = hdf.create_group('entry1') + for key, value in MOCK_METADATA.items(): + ds.create_dataset(key, data=np.array([value.encode('utf-8')])) + for key, (value, dtype) in MOCK_META_TYPED.items(): + if type(value) is list: + ds.create_dataset(key, data=np.array(value), dtype=dtype) + else: + ds.create_dataset(key, data=np.array([value]), dtype=dtype) + + ds.create_dataset('Amor/instrument_control_parameters/mu', np.array([mu]), dtype=np.float64) + ds.create_dataset('Amor/instrument_control_parameters/nu', np.array([nu]), dtype=np.float64) + + return hdf + +def compare_with_real_data(fname): + hdf = h5py.File(fname, 'r') + ds = hdf['entry1'] + for key, value in MOCK_METADATA.items(): + try: + ds[key][0].decode('utf-8') + except KeyError: + print(f'/entry1/{key} does not exist in file') + for key, (value, dtype) in MOCK_META_TYPED.items(): + try: + item = ds[key] + except KeyError: + print(f'/entry1/{key} does not exist in file') + else: + if item.dtype != dtype: + print(f'/entry1/{key} does not match {dtype}, dataset is {item.dtype}')