From d810715bc0bfcdba0ce036fa585252f606a6a23d Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 09:22:38 +0100 Subject: [PATCH 01/17] Add requirements.txt to define setup of virtual environment --- .gitignore | 2 +- requirements.txt | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) create mode 100644 requirements.txt diff --git a/.gitignore b/.gitignore index fa997c4..ae74ea7 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,4 @@ *.bak __py_cache__ raw - +.idea diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..c92ee06 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +numpy +h5py +orsopy From 7a8b11b3c79bf3b6e2bc915232486f0b5d93e41f Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 09:24:53 +0100 Subject: [PATCH 02/17] Minor format convention, no utf-8 line needed in python 3.x --- neos.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/neos.py b/neos.py index a1ddb8f..59c57d2 100644 --- a/neos.py +++ b/neos.py @@ -1,5 +1,4 @@ #!/usr/bin/env python -#-*- coding: utf-8 -*- """ eos reduces measurements performed on Amor@SINQ, PSI @@ -37,10 +36,10 @@ import platform # - format of 'call' + add '-Y' if not supplied #===================================================================================================== def commandLineArgs(): - ''' - Process command line argument. + """ + Process command line argument. The type of the default values is used for conversion and validation. - ''' + """ msg = "eos reads data from (one or several) raw file(s) of the .hdf format, \ performs various corrections, conversations and projections and exports\ the resulting reflectivity in an orso-compatible format." @@ -153,12 +152,12 @@ def commandLineArgs(): return clas.parse_args() #===================================================================================================== class Defs: - '''definition of a series of fixed parameters and constants''' + """definition of a series of fixed parameters and constants""" hdm = 6.626176e-34/1.674928e-27 # h / m lamdaCut = 2.5 # Aa #===================================================================================================== class Header: - '''orso compatible output file header content''' + """orso compatible output file header content""" def __init__(self): self.owner = None @@ -203,7 +202,7 @@ class Header: #===================================================================================================== class AmorData: - '''read meta-data and event streams from .hdf file(s), apply filters and conversions''' + """read meta-data and event streams from .hdf file(s), apply filters and conversions""" #------------------------------------------------------------------------------------------------- def __init__(self): global startTime @@ -251,7 +250,7 @@ class AmorData: return f'{path}/{fileName}' #------------------------------------------------------------------------------------------------- def expand_file_list(self, short_notation): - '''Evaluate string entry for file number lists''' + """Evaluate string entry for file number lists""" #log().debug('Executing get_flist') file_list=[] for i in short_notation.split(','): @@ -267,7 +266,7 @@ class AmorData: return sorted(file_list) #------------------------------------------------------------------------------------------------- def resolve_pixels(self): - '''determine spatial coordinats and angles from pixel number''' + """determine spatial coordinats and angles from pixel number""" det = Detector() nPixel = det.nWires * det.nStripes * det.nBlades pixelID = np.arange(nPixel) @@ -486,7 +485,7 @@ class Detector: self.distance = 4000. # mm distance from focal point to leading blade edge #===================================================================================================== def expand_file_list(short_notation): - '''Evaluate string entry for file number lists''' + """Evaluate string entry for file number lists""" #log().debug('Executing get_flist') file_list=[] for i in short_notation.split(','): From b842eaa418347cc0aac8f968bd3affd0bc068d17 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 11:13:00 +0100 Subject: [PATCH 03/17] Introduce loggers instead of print statements, options -v,-vv for verbose console output --- neos.py | 99 +++++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 72 insertions(+), 27 deletions(-) diff --git a/neos.py b/neos.py index 59c57d2..8869c86 100644 --- a/neos.py +++ b/neos.py @@ -23,6 +23,7 @@ import os import sys import subprocess import logging +import logging.handlers import argparse import numpy as np from datetime import datetime @@ -148,7 +149,11 @@ def commandLineArgs(): overwrite.add_argument("-sm", "--sampleModel", type = str, help = "1-line orso sample model description") - + + misc = clas.add_argument_group('misc') + misc.add_argument('-v', '--verbose', action='store_true') + misc.add_argument('-vv', '--debug', action='store_true') + return clas.parse_args() #===================================================================================================== class Defs: @@ -290,7 +295,7 @@ class AmorData: if self.readHeaderInfo: # read general information and first data set - print(f'# meta data from: {self.file_list[0]}') + logging.info(f'# meta data from: {self.file_list[0]}') self.hdf = h5py.File(self.file_list[0], 'r', swmr=True) title = self.hdf['entry1/title'][0].decode('utf-8') @@ -341,7 +346,7 @@ class AmorData: self.chopperDetectorDistance = self.detectorDistance - float(np.take(self.hdf['entry1/Amor/chopper/distance'], 0)) self.tofCut = defs.lamdaCut * self.chopperDetectorDistance / defs.hdm * 1.e-13 - print(f'# data from file: {fileName}') + logging.info(f'# data from file: {fileName}') try: self.mu = float(np.take(self.hdf['/entry1/Amor/master_parameters/mu/value'], 0)) self.nu = float(np.take(self.hdf['/entry1/Amor/master_parameters/nu/value'], 0)) @@ -397,7 +402,7 @@ class AmorData: header.measurement_additional_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) else: header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) - print(f'# mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kap:6.3f}') + logging.info(f'# mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kap:6.3f}') # TODO: should extract monitor from counts or beam current times time self.monitor1 = self.ctime @@ -419,7 +424,7 @@ class AmorData: if not self.startTime and not norm: self.startTime = wallTime_e[0] wallTime_e -= self.startTime - #print(f'wall time from {wallTime_e[0]} to {wallTime_e[-1]}') + logging.debug(f'wall time from {wallTime_e[0]} to {wallTime_e[-1]}') # filter 'strange' tof times > 2 tau if True: @@ -428,7 +433,7 @@ class AmorData: pixelID_e = pixelID_e[filter_e] wallTime_e = wallTime_e[filter_e] if np.shape(filter_e)[0]-np.shape(tof_e)[0] > 0.5 : - print(f'# strange times: {np.shape(filter_e)[0]-np.shape(tof_e)[0]}') + logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(tof_e)[0]}') tof_e = np.remainder( tof_e - self.tofCut + self.tau, self.tau) + self.tofCut # tof shifted to 1 frame tof_e = tof_e + self.tau * clas.chopperPhaseOffset / 180. # correction for time offset between chopper pulse and tof zero @@ -470,7 +475,7 @@ class AmorData: self.lamda_e = lamda_e[mask_e] self.wallTime_e = wallTime_e[mask_e] - print(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') + logging.info(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') #===================================================================================================== class Detector: def __init__(self): @@ -508,7 +513,7 @@ def normalisation_map(short_notation): for i in range(1, len(normalisation_list), 1): name = f'{name}_{normalisation_list[i]}' if os.path.exists(f'{clas.dataPath}/{name}.norm'): - print(f'# normalisation matrix: found and using {clas.dataPath}/{name}.norm') + logging.info(f'# normalisation matrix: found and using {clas.dataPath}/{name}.norm') norm_lz = np.loadtxt(f'{clas.dataPath}/{name}.norm') fh = open(f'{clas.dataPath}/{name}.norm', 'r') fh.readline() @@ -519,7 +524,7 @@ def normalisation_map(short_notation): normFileList[i] = entry.split('/')[-1] header.measurement_additional_files = normFileList else: - print(f'# normalisation matrix: using the files {normalisation_list}') + logging.info(f'# normalisation matrix: using the files {normalisation_list}') fromHDF.read_data(short_notation, norm=True) normAngle = fromHDF.nu - fromHDF.mu lamda_e = fromHDF.lamda_e @@ -643,7 +648,7 @@ def autoscale(q_q, R_q, dR_q, pR_q=[], pdR_q=[]): if len(filter_q[filter_q]) > 0: scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) else: - print(f'# automatic scaling not possible') + logging.warning(f'# automatic scaling not possible') scale = 1. else: filter_q = np.where(np.isnan(pR_q*R_q), False, True) @@ -653,11 +658,11 @@ def autoscale(q_q, R_q, dR_q, pR_q=[], pdR_q=[]): scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) else: - print(f'# automatic scaling not possible') + logging.warning(f'# automatic scaling not possible') scale = 1. R_q /= scale dR_q /= scale - #print(f'# scaling factor = {scale}') + logging.debug(f'# scaling factor = {scale}') return R_q, dR_q #===================================================================================================== @@ -686,7 +691,7 @@ class Grid: a, b = np.histogram([clas.qResolution], bins = resolutions) dqdq = np.matmul(b[:-1],a) if dqdq != clas.qResolution: - print(f'# changed resolution to {dqdq}') + logging.info(f'# changed resolution to {dqdq}') qq = 0.01 # linear up to qq q_grid = np.arange(0, qq, qq*dqdq) @@ -724,9 +729,47 @@ class Grid: return -np.flip(delta_grid) + 0.5*self.det.nBlades * bladeAngle #===================================================================================================== +def setup_logging(): + logger = logging.getLogger() # logging.getLogger('quicknxs') + logger.setLevel(logging.DEBUG) + # rename levels to make clear warning is can be a normal message + logging.addLevelName(logging.INFO, 'VERB') + logging.addLevelName(logging.WARNING, 'MESG') + # setting up a logger for console output + console = logging.StreamHandler(sys.__stdout__) + console.name = 'console' + formatter = logging.Formatter('# %(message)s') + console.setFormatter(formatter) + console.setLevel(logging.WARNING) + logger.addHandler(console) + + # if os.path.exists('amor_eos.log'): + # rollover = True + # else: + # rollover = False + logfile = logging.handlers.RotatingFileHandler('amor_eos.log', encoding='utf8', mode='w', + maxBytes=200*1024**2, backupCount=20) + # if rollover: logfile.doRollover() + formatter = logging.Formatter( + '[%(levelname).4s] - %(asctime)s - %(filename)s:%(lineno)i:%(funcName)s %(message)s', + '') + logfile.setFormatter(formatter) + logfile.setLevel(logging.DEBUG) + logger.addHandler(logfile) +#===================================================================================================== def main(): + setup_logging() + global startTime, grid, clas, header clas = commandLineArgs() + if clas.verbose: + logging.getLogger().handlers[0].setLevel(logging.INFO) + if clas.debug: + console = logging.getLogger().handlers[0] + console.setLevel(logging.DEBUG) + formatter = logging.Formatter('%(levelname).1s %(message)s') + console.setFormatter(formatter) + grid = Grid() header = Header() startTime = 0 @@ -734,7 +777,7 @@ def main(): os.system(f'mkdir {clas.dataPath}') fromHDF = AmorData() - print('\n######## eos - data reduction for Amor ########') + logging.warning('\n######## eos - data reduction for Amor ########') # load or create normalisation matrix if clas.normalisationFileIdentifier: @@ -746,13 +789,13 @@ def main(): norm_lz = grid.lz() normAngle = 1. - print('# normalisation matrix: none requested') + logging.warning('normalisation matrix: none requested') # load R(q_z) curve to be subtracted: if clas.subtract: sq_q, sR_q, sdR_q, sFileName = loadRqz(clas.subtract) subtract = True - print(f'# loaded background file: {sFileName}') + logging.warning(f'loaded background file: {sFileName}') header.reduction.corrections.append(f'background from \'{sFileName}\' subtracted') else: subtract = False @@ -761,7 +804,7 @@ def main(): datasetsRqz = [] datasetsRlt = [] for i, short_notation in enumerate(clas.fileIdentifier): - print('# reading input:') + logging.warning('reading input:') header.measurement_data_files = [] fromHDF.read_data(short_notation) @@ -778,9 +821,11 @@ def main(): try: stop = clas.timeSlize[2] except: - stop = wallTime_e[-1] + stop = wallTime_e[-1] + # make overwriting log lines possible by removing newline at the end + logging.StreamHandler.terminator="\r" for i, time in enumerate(np.arange(start, stop, interval)): - print(f'# time slize {i:4d}', end='\r') + logging.warning(f' time slize {i:4d}') filter_e = np.where((time < wallTime_e) & (wallTime_e < time+interval), True, False) lamda_e = fromHDF.lamda_e[filter_e] @@ -804,7 +849,7 @@ def main(): dR_q = np.sqrt( dR_q**2 + sdR_q**2 ) else: subtract = False - print(f'# background file {sFileName} not compatible with q_z scale ({len(sq_q)} vs. {len(q_q)})') + logging.warning(f'background file {sFileName} not compatible with q_z scale ({len(sq_q)} vs. {len(q_q)})') tme_q = np.ones(np.shape(q_q))*time data = np.array([q_q, R_q, dR_q, dq_q, tme_q]).T @@ -813,8 +858,9 @@ def main(): # make a copy of the header for the next iteration headerRqz = fileio.Orso.from_dict(headerRqz.to_dict()) datasetsRqz.append(orso_data) - print('') - + # reset normal logging behavior + logging.StreamHandler.terminator="\n" + logging.warning(f' time slize {i:4d}, done') else: lamda_e = fromHDF.lamda_e detZ_e = fromHDF.detZ_e @@ -853,7 +899,7 @@ def main(): R_q -= sR_q dR_q = np.sqrt( dR_q**2 + sdR_q**2 ) else: - print(f'# backgroung file {sFileName} not compatible with q_z scale ({len(sq_q)} vs. {len(q_q)})') + logging.warning(f'backgroung file {sFileName} not compatible with q_z scale ({len(sq_q)} vs. {len(q_q)})') data = np.array([q_q, R_q, dR_q, dq_q]).T orso_data = fileio.OrsoDataset(headerRqz, data) @@ -903,19 +949,18 @@ def main(): j += 1 # output - print('# output:') + logging.warning('output:') if 'Rqz.ort' in output_format_list(clas.outputFormat): - print(f'# {clas.dataPath}/{clas.outputName}.Rqz.ort') + logging.warning(f' {clas.dataPath}/{clas.outputName}.Rqz.ort') theSecondLine = f' {header.experiment.title} | {header.experiment.start_date} | sample {header.sample.name} | R(q_z)' fileio.save_orso(datasetsRqz, f'{clas.dataPath}/{clas.outputName}.Rqz.ort', data_separator='\n', comment=theSecondLine) if 'Rlt.ort' in output_format_list(clas.outputFormat): - print(f'# {clas.dataPath}/{clas.outputName}.Rlt.ort') + logging.warning(f' {clas.dataPath}/{clas.outputName}.Rlt.ort') theSecondLine = f' {header.experiment.title} | {header.experiment.start_date} | sample {header.sample.name} | R(lambda, theta)' fileio.save_orso(datasetsRlt, f'{clas.dataPath}/{clas.outputName}.Rlt.ort', data_separator='\n', comment=theSecondLine) - print('') #===================================================================================================== if __name__ == '__main__': main() From a18b82071218d453fc1434a8de2759d8c491817a Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 12:10:34 +0100 Subject: [PATCH 04/17] Separate code into multiple modules within libeos package --- .gitignore | 1 + libeos/__init__.py | 6 + libeos/command_line.py | 155 +++++++++ libeos/const.py | 6 + libeos/dataset.py | 349 +++++++++++++++++++ libeos/instrument.py | 69 ++++ libeos/logconfig.py | 42 +++ libeos/reduction.py | 158 +++++++++ neos.py | 775 +++-------------------------------------- 9 files changed, 826 insertions(+), 735 deletions(-) create mode 100644 libeos/__init__.py create mode 100644 libeos/command_line.py create mode 100644 libeos/const.py create mode 100644 libeos/dataset.py create mode 100644 libeos/instrument.py create mode 100644 libeos/logconfig.py create mode 100644 libeos/reduction.py diff --git a/.gitignore b/.gitignore index ae74ea7..0fbb9a8 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ __py_cache__ raw .idea +test_data diff --git a/libeos/__init__.py b/libeos/__init__.py new file mode 100644 index 0000000..1099066 --- /dev/null +++ b/libeos/__init__.py @@ -0,0 +1,6 @@ +""" +Package to handle data redction at AMOR instrument to be used by eos.py script. +""" + +__version__ = '2.0' +__date__ = '2024-03-04' diff --git a/libeos/command_line.py b/libeos/command_line.py new file mode 100644 index 0000000..9e639aa --- /dev/null +++ b/libeos/command_line.py @@ -0,0 +1,155 @@ +import argparse +from datetime import datetime + + +def commandLineArgs(): + """ + Process command line argument. + The type of the default values is used for conversion and validation. + """ + msg = "eos reads data from (one or several) raw file(s) of the .hdf format, \ + performs various corrections, conversations and projections and exports\ + the resulting reflectivity in an orso-compatible format." + clas = argparse.ArgumentParser(description = msg) + + input_data = clas.add_argument_group('input data') + input_data.add_argument("-n", "--fileIdentifier", + default = ['0'], + nargs = '+', + help = "file number(s) or offset (if negative)") + input_data.add_argument("-r", "--normalisationFileIdentifier", + default = [], + nargs = '+', + help = "file number(s) of normalisation measurement") + input_data.add_argument("-d", "--dataPath", + type = str, + default = '.', + help = "relative path to directory with .hdf files") + input_data.add_argument("-Y", "--year", + default = datetime.year, + type = int, + help = "year the measurement was performed") + input_data.add_argument("-sub", "--subtract", + help = "R(q_z) curve to be subtracted (in .Rqz.ort format)") + + output = clas.add_argument_group('output') + output.add_argument("-o", "--outputName", + default = "fromEOS", + help = "output file name (withot suffix)") + output.add_argument("-of", "--outputFormat", + nargs = '+', + default = ['Rqz.ort']) + output.add_argument("--offSpecular", + type = bool, + default = False, + ) + output.add_argument("-a", "--qResolution", + default = 0.01, + type = float, + help = "q_z resolution") + output.add_argument("-ts", "--timeSlize", + nargs = '+', + type = float, + help = "time slizing ,[ [,stop]]") + output.add_argument("-s", "--scale", + nargs = '+', + default = [1], + type = float, + help = "scaling factor for R(q_z)") + output.add_argument("-S", "--autoscale", + nargs = 2, + type = float, + help = "scale to 1 in the given q_z range") + + masks = clas.add_argument_group('masks') + masks.add_argument("-l", "--lambdaRange", + default = [2., 15.], + nargs = 2, + type = float, + help = "wavelength range") + masks.add_argument("-t", "--thetaRange", + default = [-12., 12.], + nargs = 2, + type = float, + help = "absolute theta range") + masks.add_argument("-T", "--thetaRangeR", + default = [-12., 12.], + nargs = 2, + type = float, + help = "relative theta range") + masks.add_argument("-y", "--yRange", + default = [11, 41], + nargs = 2, + type = int, + help = "detector y range") + masks.add_argument("-q", "--qzRange", + default = [0.005, 0.30], + nargs = 2, + type = float, + help = "q_z range") + + overwrite = clas.add_argument_group('overwrite') + overwrite.add_argument("-cs", "--chopperSpeed", + type = float, + help = "chopper speed in rpm") + overwrite.add_argument("-cp", "--chopperPhase", + default = -13.5, + type = float, + help = "chopper phase") + overwrite.add_argument("-co", "--chopperPhaseOffset", + default = -5, + type = float, + help = "phase offset between chopper opening and trigger pulse") + overwrite.add_argument("-m", "--muOffset", + default = 0., + type = float, + help = "mu offset") + overwrite.add_argument("-mu", "--mu", + default = 0, + type = float, + help ="value of mu") + overwrite.add_argument("-nu", "--nu", + default = 0, + type = float, + help = "value of nu") + overwrite.add_argument("-sm", "--sampleModel", + type = str, + help = "1-line orso sample model description") + + misc = clas.add_argument_group('misc') + misc.add_argument('-v', '--verbose', action='store_true') + misc.add_argument('-vv', '--debug', action='store_true') + + return clas.parse_args() + + +def expand_file_list(short_notation): + """Evaluate string entry for file number lists""" + #log().debug('Executing get_flist') + 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 sorted(file_list) + + +def output_format_list(outputFormat): + format_list = [] + if 'ort' in outputFormat or 'Rqz.ort' in outputFormat or 'Rqz' in outputFormat: + format_list.append('Rqz.ort') + if 'ort' in outputFormat or 'Rlt.ort' in outputFormat or 'Rlt' in outputFormat: + format_list.append('Rlt.ort') + if 'orb' in outputFormat or 'Rqz.orb' in outputFormat or 'Rqz' in outputFormat: + format_list.append('Rqz.orb') + if 'orb' in outputFormat or 'Rlt.orb' in outputFormat or 'Rlt' in outputFormat: + format_list.append('Rlt.orb') + + return sorted(format_list, reverse=True) diff --git a/libeos/const.py b/libeos/const.py new file mode 100644 index 0000000..27cd24a --- /dev/null +++ b/libeos/const.py @@ -0,0 +1,6 @@ +""" +Constants used in data reduction. +""" + +hdm = 6.626176e-34/1.674928e-27 # h / m +lamdaCut = 2.5 # Aa diff --git a/libeos/dataset.py b/libeos/dataset.py new file mode 100644 index 0000000..cf02eba --- /dev/null +++ b/libeos/dataset.py @@ -0,0 +1,349 @@ +import logging +import os +import platform +import subprocess +import sys +from datetime import datetime +from dataclasses import dataclass +from typing import Optional, Tuple + +import h5py +import numpy as np +from orsopy import fileio + +from . import __version__, const +from .instrument import Detector + +@dataclass +class DataReaderConfig: + year: int + dataPath: str + sampleModel: str + + chopperPhase: float + yRange: Tuple[float, float] + lambdaRange: Tuple[float, float] + qzRange: Tuple[float, float] + + chopperPhaseOffset: float = 0.0 + mu: Optional[float] = None + nu: Optional[float] = None + muOffset: Optional[float] = None + offSpecular: bool = False + +class Header: + """orso compatible output file header content""" + + def __init__(self): + self.owner = None + self.experiment = None + self.sample = None + self.measurement_instrument_settings = None + self.measurement_scheme = None + self.measurement_data_files = [] + self.measurement_additional_files = [] + + self.reduction = fileio.Reduction( + software = fileio.Software('eos', version=__version__), + call = ' '.join(sys.argv), + computer = platform.node(), + timestamp = datetime.now(), + creator = None, + corrections = ['histogramming in lambda and alpha_f', + 'gravity'], + ) + #------------------------------------------------------------------------------------------------- + def data_source(self): + return fileio.DataSource( + self.owner, + self.experiment, + self.sample, + fileio.Measurement( + instrument_settings = self.measurement_instrument_settings, + scheme = self.measurement_scheme, + data_files = self.measurement_data_files, + additional_files = self.measurement_additional_files, + ), + ) + #------------------------------------------------------------------------------------------------- + def columns(self): + cols = [ + fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), + fileio.Column('R', '', 'specular reflectivity'), + fileio.ErrorColumn(error_of='R', error_type='uncertainty', distribution='gaussian', value_is='sigma'), + fileio.ErrorColumn(error_of='Qz', error_type='resolution', distribution='gaussian', value_is='sigma'), + ] + return cols + + +class AmorData: + """read meta-data and event streams from .hdf file(s), apply filters and conversions""" + #------------------------------------------------------------------------------------------------- + def __init__(self, startTime, header:Header, config: DataReaderConfig): + self.startTime = startTime + self.header = header + self.config = config + #------------------------------------------------------------------------------------------------- + def read_data(self, short_notation, norm=False): + self.data_file_numbers = self.expand_file_list(short_notation) + #self.year = clas.year + self.file_list = [] + for number in self.data_file_numbers: + self.file_list.append(self.path_generator(number)) + ## read specific meta data and measurement from first file + if norm: + self.readHeaderInfo = False + else: + self.readHeaderInfo = True + + _detZ_e = [] + _lamda_e = [] + _wallTime_e = [] + for file in self.file_list: + self.read_individual_data(file, norm) + _detZ_e = np.append(_detZ_e, self.detZ_e) + _lamda_e = np.append(_lamda_e, self.lamda_e) + _wallTime_e = np.append(_wallTime_e, self.wallTime_e) + self.detZ_e = _detZ_e + self.lamda_e = _lamda_e + self.wallTime_e = _wallTime_e + + #------------------------------------------------------------------------------------------------- + def path_generator(self, number): + fileName = f'amor{self.config.year}n{number:06d}.hdf' + if os.path.exists(f'{self.config.dataPath}/{fileName}'): + path = self.config.dataPath + elif os.path.exists(fileName): + path = '.' + elif os.path.exists(f'./raw/{fileName}'): + path = './raw' + elif os.path.exists(f'../raw/{fileName}'): + path = '../raw' + elif os.path.exists(f'/afs/psi.ch/project/sinqdata/{self.config.year}/amor/{int(number/1000)}/{fileName}'): + path = '/afs/psi.ch/project/sinqdata/{self.config.year}/amor/{int(number/1000)}' + else: + sys.exit(f'# ERROR: the file {fileName} is nowhere to be found!') + return f'{path}/{fileName}' + #------------------------------------------------------------------------------------------------- + def expand_file_list(self, short_notation): + """Evaluate string entry for file number lists""" + #log().debug('Executing get_flist') + 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 sorted(file_list) + #------------------------------------------------------------------------------------------------- + def resolve_pixels(self): + """determine spatial coordinats and angles from pixel number""" + det = Detector() + nPixel = det.nWires * det.nStripes * det.nBlades + pixelID = np.arange(nPixel) + (bladeNr, bPixel) = np.divmod(pixelID, det.nWires * det.nStripes) + (bZi, detYi) = np.divmod(bPixel, det.nStripes) # z index on blade, y index on detector + detZi = bladeNr * det.nWires + bZi # z index on detector + detX = bZi * det.dX # x position in detector + # detZ = det.zero - bladeNr * det.bladeZ - bZi * det.dZ # z position on detector + bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*det.bladeZ / det.distance) ) + delta = (det.nBlades/2. - bladeNr) * bladeAngle \ + - np.rad2deg( np.arctan(bZi*det.dZ / ( det.distance + bZi * det.dX) ) ) + self.delta_z = delta[detYi==1] + return np.vstack((detYi.T, detZi.T, detX.T, delta.T)).T + #return matr + #------------------------------------------------------------------------------------------------- + def read_individual_data(self, fileName, norm=False): + pixelLookUp = self.resolve_pixels() + + self.hdf = h5py.File(fileName, 'r', swmr=True) + + if self.readHeaderInfo: + # read general information and first data set + logging.info(f'# meta data from: {self.file_list[0]}') + self.hdf = h5py.File(self.file_list[0], 'r', swmr=True) + + title = self.hdf['entry1/title'][0].decode('utf-8') + proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8') + user_name = self.hdf['entry1/user/name'][0].decode('utf-8') + user_affiliation = 'unknown' + 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') + start_date = start_time.split(' ')[0] + + if self.config.sampleModel: + model = self.config.sampleModel + else: + model = None + + # assembling orso header information + self.header.owner = fileio.Person( + name = user_name, + affiliation = user_affiliation, + contact = user_email, + ) + if user_orcid: + self.header.owner.orcid = user_orcid + self.header.experiment = fileio.Experiment( + title = title, + instrument = instrumentName, + start_date = start_date, + probe = sourceProbe, + facility = source, + proposalID = proposal_id + ) + self.header.sample = fileio.Sample( + name = sampleName, + model = model, + sample_parameters = None, + ) + self.header.measurement_scheme = 'angle- and energy-dispersive' + + 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 + + logging.info(f'# data from file: {fileName}') + try: + self.mu = float(np.take(self.hdf['/entry1/Amor/master_parameters/mu/value'], 0)) + self.nu = float(np.take(self.hdf['/entry1/Amor/master_parameters/nu/value'], 0)) + self.kap = float(np.take(self.hdf['/entry1/Amor/master_parameters/kap/value'], 0)) + self.kad = float(np.take(self.hdf['/entry1/Amor/master_parameters/kad/value'], 0)) + self.div = float(np.take(self.hdf['/entry1/Amor/master_parameters/div/value'], 0)) + self.chopperSpeed = float(np.take(self.hdf['/entry1/Amor/chopper/rotation_speed/value'], 0)) + self.chopperPhase = float(np.take(self.hdf['/entry1/Amor/chopper/phase/value'], 0)) + except(KeyError, IndexError): + logging.warning(" using parameters from nicos cache") + year_date = str(self.start_date).replace('-', '/', 1) + cachePath = '/home/amor/nicosdata/amor/cache/' + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-mu/{year_date}')).split('\t')[-1] + self.mu = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-nu/{year_date}')).split('\t')[-1] + self.nu = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kap/{year_date}')).split('\t')[-1] + self.kap = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kad/{year_date}')).split('\t')[-1] + self.kad = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-div/{year_date}')).split('\t')[-1] + self.div = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-ch1_speed/{year_date}')).split('\t')[-1] + self.chopperSpeed = float(value) + self.chopperPhase = self.config.chopperPhase + self.tau = 30. / self.chopperSpeed + + if self.config.muOffset: + self.mu += self.config.muOffset + if self.config.mu: + self.mu = self.config.mu + if self.config.nu: + self.nu = self.config.nu + + # TODO: figure out real stop time.... + self.ctime=(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-1] + - self.hdf['/entry1/Amor/detector/data/event_time_zero'][0]) / 1.e9 + fileDate = datetime.fromisoformat( self.hdf['/entry1/start_time'][0].decode('utf-8') ) + + # add header content + 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, + ) + self.header.measurement_instrument_settings.mu = fileio.Value(round(self.mu, 3), 'deg', comment='sample angle to horizon') + self.header.measurement_instrument_settings.nu = fileio.Value(round(self.nu, 3), 'deg', comment='detector angle to horizon') + if norm: + self.header.measurement_additional_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) + else: + self.header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) + logging.info(f'# mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kap:6.3f}') + + # TODO: should extract monitor from counts or beam current times time + self.monitor1 = self.ctime + self.monitor2=self.monitor1 + + # read data event streams + tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 + pixelID_e = np.array(self.hdf['/entry1/Amor/detector/data/event_id'][:], dtype=int) + totalNumber = np.shape(tof_e)[0] + + wallTime_e = np.empty(totalNumber) + dataPacket_p = np.array(self.hdf['/entry1/Amor/detector/data/event_index'][:], dtype=np.uint64) + dataPacketTime_p = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.uint64)/1e9 + #for i, index in enumerate(dataPacket_p): + # wallTime_e[index:] = dataPacketTime_p[i] + for i in range(len(dataPacket_p)-1): + wallTime_e[dataPacket_p[i]:dataPacket_p[i+1]] = dataPacketTime_p[i] + wallTime_e[dataPacket_p[-1]:] = dataPacketTime_p[-1] + if not self.startTime and not norm: + self.startTime = wallTime_e[0] + wallTime_e -= self.startTime + logging.debug(f'wall time from {wallTime_e[0]} to {wallTime_e[-1]}') + + # filter 'strange' tof times > 2 tau + if True: + filter_e = (tof_e <= 2*self.tau) + tof_e = tof_e[filter_e] + pixelID_e = pixelID_e[filter_e] + wallTime_e = wallTime_e[filter_e] + if np.shape(filter_e)[0]-np.shape(tof_e)[0] > 0.5 : + logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(tof_e)[0]}') + tof_e = np.remainder( tof_e - self.tofCut + self.tau, self.tau) + self.tofCut # tof shifted to 1 frame + tof_e = tof_e + self.tau * self.config.chopperPhaseOffset / 180. # correction for time offset between chopper pulse and tof zero + + # resolve pixel ID into y and z indicees, x position and angle + (detY_e, detZ_e, detXdist_e, delta_e) = pixelLookUp[np.int_(pixelID_e)-1,:].T + + # define mask and filter y range + mask_e = (self.config.yRange[0] <= detY_e) & (detY_e <= self.config.yRange[1]) + + # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau + # TODO: check for correctness + if not self.config.offSpecular: + tof_e -= ( delta_e / 180. ) * self.tau + + # lambda + lamda_e = 1.e13 * tof_e * const.hdm / (self.chopperDetectorDistance + detXdist_e) + self.lamdaMax = const.lamdaCut + 1.e13 * self.tau * const.hdm / (self.chopperDetectorDistance + 124.) + mask_e = np.logical_and(mask_e, (self.config.lambdaRange[0] <= lamda_e) & (lamda_e <= self.config.lambdaRange[1])) + + # alpha_f + alphaF_e = self.nu - self.mu + delta_e + + # q_z + if self.config.offSpecular: + alphaI = self.kap + self.kad + self.mu + qz_e = 2 * np.pi * ( np.sin( np.deg2rad(alphaF_e) ) + np.sin( np.deg2rad( alphaI ) ) ) / lamda_e + qx_e = 2 * np.pi * ( np.cos( np.deg2rad(alphaF_e) ) - np.cos( np.deg2rad( alphaI ) ) ) / lamda_e + self.header.measurement_scheme = 'energy-dispersive', + else: + qz_e = 4 * np.pi * np.sin( np.deg2rad(alphaF_e) ) / lamda_e + # qx_e = 0. + self.header.measurement_scheme = 'angle- and energy-dispersive' + + # filter q_z range + if self.config.qzRange[1] < 0.3 and not norm: + mask_e = np.logical_and(mask_e, (self.config.qzRange[0] <= qz_e) & (qz_e <= self.config.qzRange[1])) + + self.detZ_e = detZ_e[mask_e] + self.lamda_e = lamda_e[mask_e] + self.wallTime_e = wallTime_e[mask_e] + + logging.info(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') + diff --git a/libeos/instrument.py b/libeos/instrument.py new file mode 100644 index 0000000..7bd2c77 --- /dev/null +++ b/libeos/instrument.py @@ -0,0 +1,69 @@ +""" +Classes describing the AMOR instrument configuration used during reduction. +""" + +import logging +import numpy as np + +from . import const + +class Detector: + nBlades = 14 # number of active blades in the detector + nWires = 32 # number of wires per blade + nStripes = 64 # number of stipes per blade + angle = np.deg2rad(5.1) # deg angle of incidence of the beam on the blades (def: 5.1) + dZ = 4.0*np.sin(angle) # mm height-distance of neighboring pixels on one blade + dX = 4.0*np.cos(angle) # mm depth-distance of neighboring pixels on one blace + bladeZ = 10.455 # mm distance between detector blades + zero = 0.5*nBlades*bladeZ # mm vertical center of the detector + distance = 4000. # mm distance from focal point to leading blade edge + +class Grid: + + def __init__(self, qResolution): + self.det = Detector() + self.lamdaCut = const.lamdaCut + self.dldl = 0.005 # Delta lambda / lambda + self.qResolution = qResolution + + def q(self): + resolutions = [0.005, 0.01, 0.02, 0.025, 0.04, 0.05, 0.1, 1] + a, b = np.histogram([self.qResolution], bins = resolutions) + dqdq = np.matmul(b[:-1],a) + if dqdq != self.qResolution: + logging.info(f'# changed resolution to {dqdq}') + qq = 0.01 + # linear up to qq + q_grid = np.arange(0, qq, qq*dqdq) + # exponential from qq on + q_grid = np.append(q_grid, qq*(1.+dqdq)**np.arange(int(np.log(0.3/qq)/np.log(1+dqdq)))) + return q_grid + + def lamda(self): + lamdaMax = 16 + lamdaMin = self.lamdaCut + lamda_grid = lamdaMin*(1+self.dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+self.dldl)+1)) + return lamda_grid + + def z(self): + return np.arange(self.det.nBlades*self.det.nWires+1) + + def lz(self): + return np.ones(( np.shape(self.lamda()[:-1])[0], np.shape(self.z()[:-1])[0] )) + + def delta(self, detectorDistance): + # unused for now + bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*self.det.bladeZ / detectorDistance) ) + blade_grid = np.arctan( np.arange(33) * self.det.dZ / ( detectorDistance + np.arange(33) * self.det.dX) ) + blade_grid = np.rad2deg(blade_grid) + stepWidth = blade_grid[1] - blade_grid[0] + blade_grid = blade_grid - 0.2 * stepWidth + + delta_grid = [] + for b in np.arange(self.det.nBlades-1): + delta_grid = np.concatenate((delta_grid, blade_grid), axis=None) + blade_grid = blade_grid + bladeAngle + delta_grid = delta_grid[delta_grid0, norm_lz, np.nan) + # correct for the SM reflectivity + lamda_l = grid.lamda() + theta_z = normAngle + fromHDF.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 + 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) + norm_lz = norm_lz / Rsm_lz + if len(lamda_e) > 1e6: + head = ('normalisation matrix based on the measurements\n' + f'{fromHDF.file_list}\n' + f'nu - mu = {normAngle}\n' + f'shape= {np.shape(norm_lz)} (lambda, z)\n' + f'measured at mu = {fromHDF.mu:6.3f} deg\n' + f'N(l_lambda, z) = theta(z) / sum_i=-1..1 I(l_lambda+i, z)') + head = head.replace('../', '') + head = head.replace('./', '') + head = head.replace('raw/', '') + np.savetxt(f'{dataPath}/{name}.norm', norm_lz, header = head) + normFileList = fromHDF.file_list + return norm_lz, normAngle, normFileList + + +def project_on_lz(fromHDF, norm_lz, normAngle, lamda_e, detZ_e, grid, thetaRange, thetaRangeR, lambdaRange): + # projection on lambda-z-grid + lamda_l = grid.lamda() + theta_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z + lamda_lz = (grid.lz().T*lamda_l[:-1]).T + theta_lz = grid.lz()*theta_z + + thetaN_z = fromHDF.delta_z + normAngle + thetaN_lz = np.ones(np.shape(norm_lz))*thetaN_z + thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) + + mask_lz = np.where(np.isnan(norm_lz), False, True) + mask_lz = np.logical_and(mask_lz, np.where(np.absolute(thetaN_lz)>5e-3, True, False)) + mask_lz = np.logical_and(mask_lz, np.where(np.absolute(theta_lz)>5e-3, True, False)) + if thetaRange[1]<12: + mask_lz = np.logical_and(mask_lz, np.where(theta_lz >= thetaRange[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(theta_lz <= thetaRange[1], True, False)) + if thetaRangeR[1]<12: + t0 = fromHDF.nu - fromHDF.mu + mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 >= thetaRangeR[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 <= thetaRangeR[1], True, False)) + if lambdaRange[1]<15: + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= lambdaRange[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= lambdaRange[1], True, False)) + + # gravity correction + #theta_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) + theta_lz += np.rad2deg( np.arctan( 3.07e-10 * fromHDF.detectorDistance * lamda_lz**2 ) ) + + z_z = enumerate(theta_z) + qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz + int_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (lamda_l, grid.z())) + # cut normalisation sample horizon + int_lz = np.where(mask_lz, int_lz, np.nan) + thetaF_lz = np.where(mask_lz, theta_lz, np.nan) + + ref_lz = (int_lz * np.absolute(thetaN_lz)) / (norm_lz * np.absolute(thetaF_lz)) + err_lz = ref_lz * np.sqrt( 1/(int_lz+.1) + 1/norm_lz ) + + res_lz = np.ones((np.shape(lamda_l[:-1])[0], np.shape(theta_z)[0])) * 0.022**2 + res_lz = res_lz + (0.008/theta_lz)**2 + res_lz = qz_lz * np.sqrt(res_lz) + + return qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz + + +def project_on_qz(q_lz, R_lz, dR_lz, dq_lz, norm_lz, mask_lz, grid): + q_q = grid.q() + mask_lzf = mask_lz.flatten() + q_lzf = q_lz.flatten()[mask_lzf] + R_lzf = R_lz.flatten()[mask_lzf] + dR_lzf = dR_lz.flatten()[mask_lzf] + dq_lzf = dq_lz.flatten()[mask_lzf] + norm_lzf = norm_lz.flatten()[mask_lzf] + + N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + + R_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf * R_lzf )[0] + R_q = R_q / N_q + + dR_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dR_lzf)**2 )[0] + dR_q = np.sqrt( dR_q ) / N_q + + # TODO: different error propagations for dR and dq! + N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf**2 )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + dq_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dq_lzf)**2 )[0] + dq_q = np.sqrt( dq_q / N_q ) + + q_q = 0.5 * (q_q + np.roll(q_q, 1)) + + return q_q[1:], R_q, dR_q, dq_q + + +def autoscale(q_q, R_q, dR_q, autoscale, pR_q=[], pdR_q=[]): + if len(pR_q) == 0: + filter_q = np.where((autoscale[0]<=q_q)&(q_q<=autoscale[1]), True, False) + filter_q = np.where(dR_q>0, filter_q, False) + if len(filter_q[filter_q]) > 0: + scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) + else: + logging.warning(f'# automatic scaling not possible') + scale = 1. + else: + filter_q = np.where(np.isnan(pR_q*R_q), False, True) + filter_q = np.where(R_q>0, filter_q, False) + filter_q = np.where(pR_q>0, filter_q, False) + if len(filter_q[filter_q]) > 0: + scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ + / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) + else: + logging.warning(f'# automatic scaling not possible') + scale = 1. + R_q /= scale + dR_q /= scale + logging.debug(f'# scaling factor = {scale}') + + return R_q, dR_q diff --git a/neos.py b/neos.py index 8869c86..85fb92b 100644 --- a/neos.py +++ b/neos.py @@ -16,656 +16,28 @@ conventions (not strictly followed, yet): - to come """ -__version__ = '2.0' -__date__ = '2024-03-01' import os import sys -import subprocess import logging import logging.handlers -import argparse import numpy as np -from datetime import datetime -import h5py from orsopy import fileio -import platform + +import libeos.reduction +from libeos.command_line import commandLineArgs, output_format_list +from libeos.dataset import AmorData, DataReaderConfig, Header +from libeos.instrument import Grid +from libeos.logconfig import setup_logging, update_loglevel +from libeos.reduction import autoscale, normalisation_map, project_on_lz, project_on_qz + + #===================================================================================================== # TODO: # - calculate resolution using the chopperPhase # - deal with background correction # - format of 'call' + add '-Y' if not supplied #===================================================================================================== -def commandLineArgs(): - """ - Process command line argument. - The type of the default values is used for conversion and validation. - """ - msg = "eos reads data from (one or several) raw file(s) of the .hdf format, \ - performs various corrections, conversations and projections and exports\ - the resulting reflectivity in an orso-compatible format." - clas = argparse.ArgumentParser(description = msg) - - input_data = clas.add_argument_group('input data') - input_data.add_argument("-n", "--fileIdentifier", - default = ['0'], - nargs = '+', - help = "file number(s) or offset (if negative)") - input_data.add_argument("-r", "--normalisationFileIdentifier", - default = [], - nargs = '+', - help = "file number(s) of normalisation measurement") - input_data.add_argument("-d", "--dataPath", - type = str, - default = '.', - help = "relative path to directory with .hdf files") - input_data.add_argument("-Y", "--year", - default = datetime.year, - type = int, - help = "year the measurement was performed") - input_data.add_argument("-sub", "--subtract", - help = "R(q_z) curve to be subtracted (in .Rqz.ort format)") - - output = clas.add_argument_group('output') - output.add_argument("-o", "--outputName", - default = "fromEOS", - help = "output file name (withot suffix)") - output.add_argument("-of", "--outputFormat", - nargs = '+', - default = ['Rqz.ort']) - output.add_argument("--offSpecular", - type = bool, - default = False, - ) - output.add_argument("-a", "--qResolution", - default = 0.01, - type = float, - help = "q_z resolution") - output.add_argument("-ts", "--timeSlize", - nargs = '+', - type = float, - help = "time slizing ,[ [,stop]]") - output.add_argument("-s", "--scale", - nargs = '+', - default = [1], - type = float, - help = "scaling factor for R(q_z)") - output.add_argument("-S", "--autoscale", - nargs = 2, - type = float, - help = "scale to 1 in the given q_z range") - - masks = clas.add_argument_group('masks') - masks.add_argument("-l", "--lambdaRange", - default = [2., 15.], - nargs = 2, - type = float, - help = "wavelength range") - masks.add_argument("-t", "--thetaRange", - default = [-12., 12.], - nargs = 2, - type = float, - help = "absolute theta range") - masks.add_argument("-T", "--thetaRangeR", - default = [-12., 12.], - nargs = 2, - type = float, - help = "relative theta range") - masks.add_argument("-y", "--yRange", - default = [11, 41], - nargs = 2, - type = int, - help = "detector y range") - masks.add_argument("-q", "--qzRange", - default = [0.005, 0.30], - nargs = 2, - type = float, - help = "q_z range") - - overwrite = clas.add_argument_group('overwrite') - overwrite.add_argument("-cs", "--chopperSpeed", - type = float, - help = "chopper speed in rpm") - overwrite.add_argument("-cp", "--chopperPhase", - default = -13.5, - type = float, - help = "chopper phase") - overwrite.add_argument("-co", "--chopperPhaseOffset", - default = -5, - type = float, - help = "phase offset between chopper opening and trigger pulse") - overwrite.add_argument("-m", "--muOffset", - default = 0., - type = float, - help = "mu offset") - overwrite.add_argument("-mu", "--mu", - default = 0, - type = float, - help ="value of mu") - overwrite.add_argument("-nu", "--nu", - default = 0, - type = float, - help = "value of nu") - overwrite.add_argument("-sm", "--sampleModel", - type = str, - help = "1-line orso sample model description") - - misc = clas.add_argument_group('misc') - misc.add_argument('-v', '--verbose', action='store_true') - misc.add_argument('-vv', '--debug', action='store_true') - - return clas.parse_args() -#===================================================================================================== -class Defs: - """definition of a series of fixed parameters and constants""" - hdm = 6.626176e-34/1.674928e-27 # h / m - lamdaCut = 2.5 # Aa -#===================================================================================================== -class Header: - """orso compatible output file header content""" - - def __init__(self): - self.owner = None - self.experiment = None - self.sample = None - self.measurement_instrument_settings = None - self.measurement_scheme = None - self.measurement_data_files = [] - self.measurement_additional_files = [] - - self.reduction = fileio.Reduction( - software = fileio.Software('eos', version=__version__), - call = ' '.join(sys.argv), - computer = platform.node(), - timestamp = datetime.now(), - creator = None, - corrections = ['histogramming in lambda and alpha_f', - 'gravity'], - ) - #------------------------------------------------------------------------------------------------- - def data_source(self): - return fileio.DataSource( - self.owner, - self.experiment, - self.sample, - fileio.Measurement( - instrument_settings = self.measurement_instrument_settings, - scheme = self.measurement_scheme, - data_files = self.measurement_data_files, - additional_files = self.measurement_additional_files, - ), - ) - #------------------------------------------------------------------------------------------------- - def columns(self): - cols = [ - fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), - fileio.Column('R', '', 'specular reflectivity'), - fileio.ErrorColumn(error_of='R', error_type='uncertainty', distribution='gaussian', value_is='sigma'), - fileio.ErrorColumn(error_of='Qz', error_type='resolution', distribution='gaussian', value_is='sigma'), - ] - return cols - -#===================================================================================================== -class AmorData: - """read meta-data and event streams from .hdf file(s), apply filters and conversions""" - #------------------------------------------------------------------------------------------------- - def __init__(self): - global startTime - self.startTime = startTime - #------------------------------------------------------------------------------------------------- - def read_data(self, short_notation, norm=False): - self.data_file_numbers = self.expand_file_list(short_notation) - #self.year = clas.year - self.file_list = [] - for number in self.data_file_numbers: - self.file_list.append(self.path_generator(number)) - ## read specific meta data and measurement from first file - if norm: - self.readHeaderInfo = False - else: - self.readHeaderInfo = True - - _detZ_e = [] - _lamda_e = [] - _wallTime_e = [] - for file in self.file_list: - self.read_individual_data(file, norm) - _detZ_e = np.append(_detZ_e, self.detZ_e) - _lamda_e = np.append(_lamda_e, self.lamda_e) - _wallTime_e = np.append(_wallTime_e, self.wallTime_e) - self.detZ_e = _detZ_e - self.lamda_e = _lamda_e - self.wallTime_e = _wallTime_e - - #------------------------------------------------------------------------------------------------- - def path_generator(self, number): - fileName = f'amor{clas.year}n{number:06d}.hdf' - if os.path.exists(f'{clas.dataPath}/{fileName}'): - path = clas.dataPath - elif os.path.exists(fileName): - path = '.' - elif os.path.exists(f'./raw/{fileName}'): - path = './raw' - elif os.path.exists(f'../raw/{fileName}'): - path = '../raw' - elif os.path.exists(f'/afs/psi.ch/project/sinqdata/{clas.year}/amor/{int(number/1000)}/{fileName}'): - path = '/afs/psi.ch/project/sinqdata/{clas.year}/amor/{int(number/1000)}' - else: - sys.exit(f'# ERROR: the file {fileName} is nowhere to be found!') - return f'{path}/{fileName}' - #------------------------------------------------------------------------------------------------- - def expand_file_list(self, short_notation): - """Evaluate string entry for file number lists""" - #log().debug('Executing get_flist') - 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 sorted(file_list) - #------------------------------------------------------------------------------------------------- - def resolve_pixels(self): - """determine spatial coordinats and angles from pixel number""" - det = Detector() - nPixel = det.nWires * det.nStripes * det.nBlades - pixelID = np.arange(nPixel) - (bladeNr, bPixel) = np.divmod(pixelID, det.nWires * det.nStripes) - (bZi, detYi) = np.divmod(bPixel, det.nStripes) # z index on blade, y index on detector - detZi = bladeNr * det.nWires + bZi # z index on detector - detX = bZi * det.dX # x position in detector - # detZ = det.zero - bladeNr * det.bladeZ - bZi * det.dZ # z position on detector - bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*det.bladeZ / det.distance) ) - delta = (det.nBlades/2. - bladeNr) * bladeAngle \ - - np.rad2deg( np.arctan(bZi*det.dZ / ( det.distance + bZi * det.dX) ) ) - self.delta_z = delta[detYi==1] - return np.vstack((detYi.T, detZi.T, detX.T, delta.T)).T - #return matr - #------------------------------------------------------------------------------------------------- - def read_individual_data(self, fileName, norm=False): - defs = Defs() - pixelLookUp = self.resolve_pixels() - - self.hdf = h5py.File(fileName, 'r', swmr=True) - - if self.readHeaderInfo: - # read general information and first data set - logging.info(f'# meta data from: {self.file_list[0]}') - self.hdf = h5py.File(self.file_list[0], 'r', swmr=True) - - title = self.hdf['entry1/title'][0].decode('utf-8') - proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8') - user_name = self.hdf['entry1/user/name'][0].decode('utf-8') - user_affiliation = 'unknown' - 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') - start_date = start_time.split(' ')[0] - - if clas.sampleModel: - model = clas.sampleModel - else: - model = None - - # assembling orso header information - header.owner = fileio.Person( - name = user_name, - affiliation = user_affiliation, - contact = user_email, - ) - if user_orcid: - header.owner.orcid = user_orcid - header.experiment = fileio.Experiment( - title = title, - instrument = instrumentName, - start_date = start_date, - probe = sourceProbe, - facility = source, - proposalID = proposal_id - ) - header.sample = fileio.Sample( - name = sampleName, - model = model, - sample_parameters = None, - ) - header.measurement_scheme = 'angle- and energy-dispersive' - - 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 = defs.lamdaCut * self.chopperDetectorDistance / defs.hdm * 1.e-13 - - logging.info(f'# data from file: {fileName}') - try: - self.mu = float(np.take(self.hdf['/entry1/Amor/master_parameters/mu/value'], 0)) - self.nu = float(np.take(self.hdf['/entry1/Amor/master_parameters/nu/value'], 0)) - self.kap = float(np.take(self.hdf['/entry1/Amor/master_parameters/kap/value'], 0)) - self.kad = float(np.take(self.hdf['/entry1/Amor/master_parameters/kad/value'], 0)) - self.div = float(np.take(self.hdf['/entry1/Amor/master_parameters/div/value'], 0)) - self.chopperSpeed = float(np.take(self.hdf['/entry1/Amor/chopper/rotation_speed/value'], 0)) - self.chopperPhase = float(np.take(self.hdf['/entry1/Amor/chopper/phase/value'], 0)) - except(KeyError, IndexError): - logging.warning(" using parameters from nicos cache") - year_date = str(self.start_date).replace('-', '/', 1) - cachePath = '/home/amor/nicosdata/amor/cache/' - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-mu/{year_date}')).split('\t')[-1] - self.mu = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-nu/{year_date}')).split('\t')[-1] - self.nu = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kap/{year_date}')).split('\t')[-1] - self.kap = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kad/{year_date}')).split('\t')[-1] - self.kad = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-div/{year_date}')).split('\t')[-1] - self.div = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-ch1_speed/{year_date}')).split('\t')[-1] - self.chopperSpeed = float(value) - self.chopperPhase = clas.chopperPhase - self.tau = 30. / self.chopperSpeed - - if clas.muOffset: - self.mu += clas.muOffset - if clas.mu: - self.mu = clas.mu - if clas.nu: - self.mu = clas.nu - - # TODO: figure out real stop time.... - self.ctime=(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-1] - - self.hdf['/entry1/Amor/detector/data/event_time_zero'][0]) / 1.e9 - fileDate = datetime.fromisoformat( self.hdf['/entry1/start_time'][0].decode('utf-8') ) - - # add header content - if self.readHeaderInfo: - self.readHeaderInfo = False - 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(defs.lamdaCut, clas.lambdaRange[1], 'angstrom'), - polarization = fileio.Polarization.unpolarized, - ) - header.measurement_instrument_settings.mu = fileio.Value(round(self.mu, 3), 'deg', comment='sample angle to horizon') - header.measurement_instrument_settings.nu = fileio.Value(round(self.nu, 3), 'deg', comment='detector angle to horizon') - if norm: - header.measurement_additional_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) - else: - header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) - logging.info(f'# mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kap:6.3f}') - - # TODO: should extract monitor from counts or beam current times time - self.monitor1 = self.ctime - self.monitor2=self.monitor1 - - # read data event streams - tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 - pixelID_e = np.array(self.hdf['/entry1/Amor/detector/data/event_id'][:], dtype=int) - totalNumber = np.shape(tof_e)[0] - - wallTime_e = np.empty(totalNumber) - dataPacket_p = np.array(self.hdf['/entry1/Amor/detector/data/event_index'][:], dtype=np.uint64) - dataPacketTime_p = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.uint64)/1e9 - #for i, index in enumerate(dataPacket_p): - # wallTime_e[index:] = dataPacketTime_p[i] - for i in range(len(dataPacket_p)-1): - wallTime_e[dataPacket_p[i]:dataPacket_p[i+1]] = dataPacketTime_p[i] - wallTime_e[dataPacket_p[-1]:] = dataPacketTime_p[-1] - if not self.startTime and not norm: - self.startTime = wallTime_e[0] - wallTime_e -= self.startTime - logging.debug(f'wall time from {wallTime_e[0]} to {wallTime_e[-1]}') - - # filter 'strange' tof times > 2 tau - if True: - filter_e = (tof_e <= 2*self.tau) - tof_e = tof_e[filter_e] - pixelID_e = pixelID_e[filter_e] - wallTime_e = wallTime_e[filter_e] - if np.shape(filter_e)[0]-np.shape(tof_e)[0] > 0.5 : - logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(tof_e)[0]}') - tof_e = np.remainder( tof_e - self.tofCut + self.tau, self.tau) + self.tofCut # tof shifted to 1 frame - tof_e = tof_e + self.tau * clas.chopperPhaseOffset / 180. # correction for time offset between chopper pulse and tof zero - - # resolve pixel ID into y and z indicees, x position and angle - (detY_e, detZ_e, detXdist_e, delta_e) = pixelLookUp[np.int_(pixelID_e)-1,:].T - - # define mask and filter y range - mask_e = (clas.yRange[0] <= detY_e) & (detY_e <= clas.yRange[1]) - - # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau - # TODO: check for correctness - if not clas.offSpecular: - tof_e -= ( delta_e / 180. ) * self.tau - - # lambda - lamda_e = 1.e13 * tof_e * defs.hdm / (self.chopperDetectorDistance + detXdist_e) - self.lamdaMax = defs.lamdaCut + 1.e13 * self.tau * defs.hdm / (self.chopperDetectorDistance + 124.) - mask_e = np.logical_and(mask_e, (clas.lambdaRange[0] <= lamda_e) & (lamda_e <= clas.lambdaRange[1])) - - # alpha_f - alphaF_e = self.nu - self.mu + delta_e - - # q_z - if clas.offSpecular: - alphaI = self.kap + self.kad + self.mu - qz_e = 2 * np.pi * ( np.sin( np.deg2rad(alphaF_e) ) + np.sin( np.deg2rad( alphaI ) ) ) / lamda_e - qx_e = 2 * np.pi * ( np.cos( np.deg2rad(alphaF_e) ) - np.cos( np.deg2rad( alphaI ) ) ) / lamda_e - header.measurement_scheme = 'energy-dispersive', - else: - qz_e = 4 * np.pi * np.sin( np.deg2rad(alphaF_e) ) / lamda_e - # qx_e = 0. - header.measurement_scheme = 'angle- and energy-dispersive' - - # filter q_z range - if clas.qzRange[1] < 0.3 and not norm: - mask_e = np.logical_and(mask_e, (clas.qzRange[0] <= qz_e) & (qz_e <= clas.qzRange[1])) - - self.detZ_e = detZ_e[mask_e] - self.lamda_e = lamda_e[mask_e] - self.wallTime_e = wallTime_e[mask_e] - - logging.info(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') -#===================================================================================================== -class Detector: - def __init__(self): - self.nBlades = 14 # number of active blades in the detector - self.nWires = 32 # number of wires per blade - self.nStripes = 64 # number of stipes per blade - self.angle = np.deg2rad(5.1) # deg angle of incidence of the beam on the blades (def: 5.1) - self.dZ = 4.0 * np.sin(self.angle) # mm height-distance of neighboring pixels on one blade - self.dX = 4.0 * np.cos(self.angle) # mm depth-distance of neighboring pixels on one blace - self.bladeZ = 10.455 # mm distance between detector blades - self.zero = 0.5 * self.nBlades * self.bladeZ # mm vertical center of the detector - self.distance = 4000. # mm distance from focal point to leading blade edge -#===================================================================================================== -def expand_file_list(short_notation): - """Evaluate string entry for file number lists""" - #log().debug('Executing get_flist') - 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 sorted(file_list) -#===================================================================================================== -def normalisation_map(short_notation): - fromHDF = AmorData() - normalisation_list = expand_file_list(short_notation) - name = str(normalisation_list[0]) - for i in range(1, len(normalisation_list), 1): - name = f'{name}_{normalisation_list[i]}' - if os.path.exists(f'{clas.dataPath}/{name}.norm'): - logging.info(f'# normalisation matrix: found and using {clas.dataPath}/{name}.norm') - norm_lz = np.loadtxt(f'{clas.dataPath}/{name}.norm') - fh = open(f'{clas.dataPath}/{name}.norm', 'r') - fh.readline() - normFileList = fh.readline().split('[')[1].split(']')[0].replace('\'', '').split(', ') - normAngle = float(fh.readline().split('= ')[1]) - fh.close() - for i, entry in enumerate(normFileList): - normFileList[i] = entry.split('/')[-1] - header.measurement_additional_files = normFileList - else: - logging.info(f'# normalisation matrix: using the files {normalisation_list}') - fromHDF.read_data(short_notation, norm=True) - normAngle = fromHDF.nu - fromHDF.mu - lamda_e = fromHDF.lamda_e - detZ_e = fromHDF.detZ_e - norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (grid.lamda(), grid.z())) - norm_lz = np.where(norm_lz>0, norm_lz, np.nan) - # correct for the SM reflectivity - lamda_l = grid.lamda() - theta_z = normAngle + fromHDF.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 - 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) - norm_lz = norm_lz / Rsm_lz - if len(lamda_e) > 1e6: - head = ('normalisation matrix based on the measurements\n' - f'{fromHDF.file_list}\n' - f'nu - mu = {normAngle}\n' - f'shape= {np.shape(norm_lz)} (lambda, z)\n' - f'measured at mu = {fromHDF.mu:6.3f} deg\n' - f'N(l_lambda, z) = theta(z) / sum_i=-1..1 I(l_lambda+i, z)') - head = head.replace('../', '') - head = head.replace('./', '') - head = head.replace('raw/', '') - np.savetxt(f'{clas.dataPath}/{name}.norm', norm_lz, header = head) - normFileList = fromHDF.file_list - return norm_lz, normAngle, normFileList -#===================================================================================================== -def output_format_list(outputFormat): - format_list = [] - if 'ort' in outputFormat or 'Rqz.ort' in outputFormat or 'Rqz' in outputFormat: - format_list.append('Rqz.ort') - if 'ort' in outputFormat or 'Rlt.ort' in outputFormat or 'Rlt' in outputFormat: - format_list.append('Rlt.ort') - if 'orb' in outputFormat or 'Rqz.orb' in outputFormat or 'Rqz' in outputFormat: - format_list.append('Rqz.orb') - if 'orb' in outputFormat or 'Rlt.orb' in outputFormat or 'Rlt' in outputFormat: - format_list.append('Rlt.orb') - - return sorted(format_list, reverse=True) -#===================================================================================================== -def project_on_lz(fromHDF, norm_lz, normAngle, lamda_e, detZ_e): - # projection on lambda-z-grid - lamda_l = grid.lamda() - theta_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z - lamda_lz = (grid.lz().T*lamda_l[:-1]).T - theta_lz = grid.lz()*theta_z - - thetaN_z = fromHDF.delta_z + normAngle - thetaN_lz = np.ones(np.shape(norm_lz))*thetaN_z - thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) - - mask_lz = np.where(np.isnan(norm_lz), False, True) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(thetaN_lz)>5e-3, True, False)) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(theta_lz)>5e-3, True, False)) - if clas.thetaRange[1]<12: - mask_lz = np.logical_and(mask_lz, np.where(theta_lz >= clas.thetaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(theta_lz <= clas.thetaRange[1], True, False)) - if clas.thetaRangeR[1]<12: - t0 = fromHDF.nu - fromHDF.mu - mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 >= clas.thetaRangeR[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 <= clas.thetaRangeR[1], True, False)) - if clas.lambdaRange[1]<15: - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= clas.lambdaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= clas.lambdaRange[1], True, False)) - - # gravity correction - #theta_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) - theta_lz += np.rad2deg( np.arctan( 3.07e-10 * fromHDF.detectorDistance * lamda_lz**2 ) ) - - z_z = enumerate(theta_z) - qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz - int_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (lamda_l, grid.z())) - # cut normalisation sample horizon - int_lz = np.where(mask_lz, int_lz, np.nan) - thetaF_lz = np.where(mask_lz, theta_lz, np.nan) - - ref_lz = (int_lz * np.absolute(thetaN_lz)) / (norm_lz * np.absolute(thetaF_lz)) - err_lz = ref_lz * np.sqrt( 1/(int_lz+.1) + 1/norm_lz ) - - res_lz = np.ones((np.shape(lamda_l[:-1])[0], np.shape(theta_z)[0])) * 0.022**2 - res_lz = res_lz + (0.008/theta_lz)**2 - res_lz = qz_lz * np.sqrt(res_lz) - - return qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz -#===================================================================================================== -def project_on_qz(q_lz, R_lz, dR_lz, dq_lz, norm_lz, mask_lz): - q_q = grid.q() - mask_lzf = mask_lz.flatten() - q_lzf = q_lz.flatten()[mask_lzf] - R_lzf = R_lz.flatten()[mask_lzf] - dR_lzf = dR_lz.flatten()[mask_lzf] - dq_lzf = dq_lz.flatten()[mask_lzf] - norm_lzf = norm_lz.flatten()[mask_lzf] - - N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - - R_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf * R_lzf )[0] - R_q = R_q / N_q - - dR_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dR_lzf)**2 )[0] - dR_q = np.sqrt( dR_q ) / N_q - - # TODO: different error propagations for dR and dq! - N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf**2 )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - dq_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dq_lzf)**2 )[0] - dq_q = np.sqrt( dq_q / N_q ) - - q_q = 0.5 * (q_q + np.roll(q_q, 1)) - - return q_q[1:], R_q, dR_q, dq_q -#===================================================================================================== -def autoscale(q_q, R_q, dR_q, pR_q=[], pdR_q=[]): - if len(pR_q) == 0: - filter_q = np.where((clas.autoscale[0]<=q_q)&(q_q<=clas.autoscale[1]), True, False) - filter_q = np.where(dR_q>0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) - else: - logging.warning(f'# automatic scaling not possible') - scale = 1. - else: - filter_q = np.where(np.isnan(pR_q*R_q), False, True) - filter_q = np.where(R_q>0, filter_q, False) - filter_q = np.where(pR_q>0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ - / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) - else: - logging.warning(f'# automatic scaling not possible') - scale = 1. - R_q /= scale - dR_q /= scale - logging.debug(f'# scaling factor = {scale}') - - return R_q, dR_q -#===================================================================================================== def loadRqz(name): if os.path.exists(f'{clas.dataPath}/{name}'): @@ -678,111 +50,40 @@ def loadRqz(name): q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) return q_q, Sq_q, dS_q, fileName -#===================================================================================================== -class Grid: - def __init__(self): - self.det = Detector() - self.lamdaCut = Defs.lamdaCut - self.dldl = 0.005 # Delta lambda / lambda - - def q(self): - resolutions = [0.005, 0.01, 0.02, 0.025, 0.04, 0.05, 0.1, 1] - a, b = np.histogram([clas.qResolution], bins = resolutions) - dqdq = np.matmul(b[:-1],a) - if dqdq != clas.qResolution: - logging.info(f'# changed resolution to {dqdq}') - qq = 0.01 - # linear up to qq - q_grid = np.arange(0, qq, qq*dqdq) - # exponential from qq on - q_grid = np.append(q_grid, qq*(1.+dqdq)**np.arange(int(np.log(0.3/qq)/np.log(1+dqdq)))) - return q_grid - - def lamda(self): - lamdaMax = 16 - lamdaMin = self.lamdaCut - lamda_grid = lamdaMin*(1+self.dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+self.dldl)+1)) - return lamda_grid - - def z(self): - return np.arange(self.det.nBlades*self.det.nWires+1) - - def lz(self): - return np.ones(( np.shape(self.lamda()[:-1])[0], np.shape(self.z()[:-1])[0] )) - - def delta(self, detectorDistance): - # unused for now - bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*self.det.bladeZ / detectorDistance) ) - blade_grid = np.arctan( np.arange(33) * self.det.dZ / ( detectorDistance + np.arange(33) * self.det.dX) ) - blade_grid = np.rad2deg(blade_grid) - stepWidth = blade_grid[1] - blade_grid[0] - blade_grid = blade_grid - 0.2 * stepWidth - - delta_grid = [] - for b in np.arange(self.det.nBlades-1): - delta_grid = np.concatenate((delta_grid, blade_grid), axis=None) - blade_grid = blade_grid + bladeAngle - delta_grid = delta_grid[delta_grid Date: Mon, 4 Mar 2024 13:31:23 +0100 Subject: [PATCH 05/17] Move all reduction code to AmorReducer class and options to dataclasses, all command line handling separated --- libeos/command_line.py | 40 +++- libeos/dataset.py | 19 +- libeos/logconfig.py | 1 + libeos/options.py | 43 ++++ libeos/reduction.py | 520 +++++++++++++++++++++++++++++------------ neos.py | 247 +------------------- 6 files changed, 465 insertions(+), 405 deletions(-) create mode 100644 libeos/options.py diff --git a/libeos/command_line.py b/libeos/command_line.py index 9e639aa..d43053d 100644 --- a/libeos/command_line.py +++ b/libeos/command_line.py @@ -1,6 +1,9 @@ import argparse from datetime import datetime +from .logconfig import update_loglevel +from .options import DataReaderConfig, OutputConfig, ReductionConfig + def commandLineArgs(): """ @@ -151,5 +154,40 @@ def output_format_list(outputFormat): format_list.append('Rqz.orb') if 'orb' in outputFormat or 'Rlt.orb' in outputFormat or 'Rlt' in outputFormat: format_list.append('Rlt.orb') - return sorted(format_list, reverse=True) + +def command_line_options(): + clas = commandLineArgs() + update_loglevel(clas.verbose, clas.debug) + + reader_config = DataReaderConfig( + year=clas.year, + dataPath=clas.dataPath, + sampleModel=clas.sampleModel, + chopperPhase=clas.chopperPhase, + chopperPhaseOffset=clas.chopperPhaseOffset, + yRange=clas.yRange, + lambdaRange=clas.lambdaRange, + qzRange=clas.qzRange, + offSpecular=clas.offSpecular, + mu=clas.mu, + nu=clas.nu, + muOffset=clas.muOffset + ) + reduction_config = ReductionConfig( + qResolution=clas.qResolution, + autoscale=clas.autoscale, + thetaRange=clas.thetaRange, + thetaRangeR=clas.thetaRangeR, + lambdaRange=clas.lambdaRange, + fileIdentifier=clas.fileIdentifier, + scale=clas.scale, + subtract=clas.subtract, + normalisationFileIdentifier=clas.normalisationFileIdentifier, + timeSlize=clas.timeSlize + ) + output_config = OutputConfig( + outputFormats=output_format_list(clas.outputFormat), + outputName=clas.outputName + ) + return reader_config, reduction_config, output_config \ No newline at end of file diff --git a/libeos/dataset.py b/libeos/dataset.py index cf02eba..10cf34d 100644 --- a/libeos/dataset.py +++ b/libeos/dataset.py @@ -4,8 +4,6 @@ import platform import subprocess import sys from datetime import datetime -from dataclasses import dataclass -from typing import Optional, Tuple import h5py import numpy as np @@ -13,23 +11,8 @@ from orsopy import fileio from . import __version__, const from .instrument import Detector +from .options import DataReaderConfig -@dataclass -class DataReaderConfig: - year: int - dataPath: str - sampleModel: str - - chopperPhase: float - yRange: Tuple[float, float] - lambdaRange: Tuple[float, float] - qzRange: Tuple[float, float] - - chopperPhaseOffset: float = 0.0 - mu: Optional[float] = None - nu: Optional[float] = None - muOffset: Optional[float] = None - offSpecular: bool = False class Header: """orso compatible output file header content""" diff --git a/libeos/logconfig.py b/libeos/logconfig.py index 031393f..6f80945 100644 --- a/libeos/logconfig.py +++ b/libeos/logconfig.py @@ -3,6 +3,7 @@ Setup for the logging of eos. """ import sys import logging +import logging.handlers def setup_logging(): logger = logging.getLogger() # logging.getLogger('quicknxs') diff --git a/libeos/options.py b/libeos/options.py new file mode 100644 index 0000000..606d6da --- /dev/null +++ b/libeos/options.py @@ -0,0 +1,43 @@ +""" +Classes for stroing various configurations needed for reduction. +""" +from dataclasses import dataclass, field +from typing import Optional, Tuple + + +@dataclass +class DataReaderConfig: + year: int + dataPath: str + sampleModel: str + + chopperPhase: float + yRange: Tuple[float, float] + lambdaRange: Tuple[float, float] + qzRange: Tuple[float, float] + + chopperPhaseOffset: float = 0.0 + mu: Optional[float] = None + nu: Optional[float] = None + muOffset: Optional[float] = None + offSpecular: bool = False + +@dataclass +class ReductionConfig: + qResolution: float + autoscale: Tuple[bool, bool] + thetaRange: Tuple[float, float] + thetaRangeR: Tuple[float, float] + lambdaRange: Tuple[float, float] + + fileIdentifier: list = field(default_factory=lambda: ["0"]) + scale: list = field(default_factory=lambda: [1]) #per file scaling, if less elements then files use the last one + + subtract: Optional[str] = None + normalisationFileIdentifier: Optional[list] = None + timeSlize: Optional[list] = None + +@dataclass +class OutputConfig: + outputFormats: list #output_format_list(clas.outputFormat) + outputName: str diff --git a/libeos/reduction.py b/libeos/reduction.py index ba1b83d..e41a526 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -1,158 +1,382 @@ import logging import os +import sys import numpy as np +from orsopy import fileio -from libeos.command_line import expand_file_list -from libeos.dataset import AmorData +from .command_line import expand_file_list +from .dataset import AmorData, Header +from .options import DataReaderConfig, ReductionConfig, OutputConfig +from .instrument import Grid +class AmorReduction: + def __init__(self, + data_reader_config: DataReaderConfig, + reduction_config: ReductionConfig, + output_config: OutputConfig): + self.data_reader_config = data_reader_config + self.reduction_config = reduction_config + self.output_config = output_config -def normalisation_map(short_notation, header, grid, dataPath): - fromHDF = AmorData() - normalisation_list = expand_file_list(short_notation) - name = str(normalisation_list[0]) - for i in range(1, len(normalisation_list), 1): - name = f'{name}_{normalisation_list[i]}' - if os.path.exists(f'{dataPath}/{name}.norm'): - logging.info(f'# normalisation matrix: found and using {dataPath}/{name}.norm') - norm_lz = np.loadtxt(f'{dataPath}/{name}.norm') - fh = open(f'{dataPath}/{name}.norm', 'r') - fh.readline() - normFileList = fh.readline().split('[')[1].split(']')[0].replace('\'', '').split(', ') - normAngle = float(fh.readline().split('= ')[1]) - fh.close() - for i, entry in enumerate(normFileList): - normFileList[i] = entry.split('/')[-1] - header.measurement_additional_files = normFileList - else: - logging.info(f'# normalisation matrix: using the files {normalisation_list}') - fromHDF.read_data(short_notation, norm=True) - normAngle = fromHDF.nu - fromHDF.mu - lamda_e = fromHDF.lamda_e - detZ_e = fromHDF.detZ_e - norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (grid.lamda(), grid.z())) - norm_lz = np.where(norm_lz>0, norm_lz, np.nan) - # correct for the SM reflectivity - lamda_l = grid.lamda() - theta_z = normAngle + fromHDF.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 - 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) - norm_lz = norm_lz / Rsm_lz - if len(lamda_e) > 1e6: - head = ('normalisation matrix based on the measurements\n' - f'{fromHDF.file_list}\n' - f'nu - mu = {normAngle}\n' - f'shape= {np.shape(norm_lz)} (lambda, z)\n' - f'measured at mu = {fromHDF.mu:6.3f} deg\n' - f'N(l_lambda, z) = theta(z) / sum_i=-1..1 I(l_lambda+i, z)') - head = head.replace('../', '') - head = head.replace('./', '') - head = head.replace('raw/', '') - np.savetxt(f'{dataPath}/{name}.norm', norm_lz, header = head) - normFileList = fromHDF.file_list - return norm_lz, normAngle, normFileList + def reduce(self): + self.grid = Grid(self.reduction_config.qResolution) + self.header = Header() + self.startTime = 0 + if not os.path.exists(f'{self.data_reader_config.dataPath}'): + os.system(f'mkdir {self.data_reader_config.dataPath}') + fromHDF = AmorData(self.startTime, header=self.header, config=self.data_reader_config) + logging.warning('\n######## eos - data reduction for Amor ########') - -def project_on_lz(fromHDF, norm_lz, normAngle, lamda_e, detZ_e, grid, thetaRange, thetaRangeR, lambdaRange): - # projection on lambda-z-grid - lamda_l = grid.lamda() - theta_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z - lamda_lz = (grid.lz().T*lamda_l[:-1]).T - theta_lz = grid.lz()*theta_z - - thetaN_z = fromHDF.delta_z + normAngle - thetaN_lz = np.ones(np.shape(norm_lz))*thetaN_z - thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) - - mask_lz = np.where(np.isnan(norm_lz), False, True) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(thetaN_lz)>5e-3, True, False)) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(theta_lz)>5e-3, True, False)) - if thetaRange[1]<12: - mask_lz = np.logical_and(mask_lz, np.where(theta_lz >= thetaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(theta_lz <= thetaRange[1], True, False)) - if thetaRangeR[1]<12: - t0 = fromHDF.nu - fromHDF.mu - mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 >= thetaRangeR[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 <= thetaRangeR[1], True, False)) - if lambdaRange[1]<15: - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= lambdaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= lambdaRange[1], True, False)) - - # gravity correction - #theta_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) - theta_lz += np.rad2deg( np.arctan( 3.07e-10 * fromHDF.detectorDistance * lamda_lz**2 ) ) - - z_z = enumerate(theta_z) - qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz - int_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (lamda_l, grid.z())) - # cut normalisation sample horizon - int_lz = np.where(mask_lz, int_lz, np.nan) - thetaF_lz = np.where(mask_lz, theta_lz, np.nan) - - ref_lz = (int_lz * np.absolute(thetaN_lz)) / (norm_lz * np.absolute(thetaF_lz)) - err_lz = ref_lz * np.sqrt( 1/(int_lz+.1) + 1/norm_lz ) - - res_lz = np.ones((np.shape(lamda_l[:-1])[0], np.shape(theta_z)[0])) * 0.022**2 - res_lz = res_lz + (0.008/theta_lz)**2 - res_lz = qz_lz * np.sqrt(res_lz) - - return qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz - - -def project_on_qz(q_lz, R_lz, dR_lz, dq_lz, norm_lz, mask_lz, grid): - q_q = grid.q() - mask_lzf = mask_lz.flatten() - q_lzf = q_lz.flatten()[mask_lzf] - R_lzf = R_lz.flatten()[mask_lzf] - dR_lzf = dR_lz.flatten()[mask_lzf] - dq_lzf = dq_lz.flatten()[mask_lzf] - norm_lzf = norm_lz.flatten()[mask_lzf] - - N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - - R_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf * R_lzf )[0] - R_q = R_q / N_q - - dR_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dR_lzf)**2 )[0] - dR_q = np.sqrt( dR_q ) / N_q - - # TODO: different error propagations for dR and dq! - N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf**2 )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - dq_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dq_lzf)**2 )[0] - dq_q = np.sqrt( dq_q / N_q ) - - q_q = 0.5 * (q_q + np.roll(q_q, 1)) - - return q_q[1:], R_q, dR_q, dq_q - - -def autoscale(q_q, R_q, dR_q, autoscale, pR_q=[], pdR_q=[]): - if len(pR_q) == 0: - filter_q = np.where((autoscale[0]<=q_q)&(q_q<=autoscale[1]), True, False) - filter_q = np.where(dR_q>0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) + # load or create normalisation matrix + if self.reduction_config.normalisationFileIdentifier: + normalise = True + norm_lz, normAngle, normFileList = self.normalisation_map(self.reduction_config.normalisationFileIdentifier[0]) + self.header.reduction.corrections.append('normalisation with \'additional files\'') else: - logging.warning(f'# automatic scaling not possible') - scale = 1. - else: - filter_q = np.where(np.isnan(pR_q*R_q), False, True) - filter_q = np.where(R_q>0, filter_q, False) - filter_q = np.where(pR_q>0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ - / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) - else: - logging.warning(f'# automatic scaling not possible') - scale = 1. - R_q /= scale - dR_q /= scale - logging.debug(f'# scaling factor = {scale}') + normalise = False + norm_lz = self.grid.lz() + normAngle = 1. + + logging.warning('normalisation matrix: none requested') + + # load R(q_z) curve to be subtracted: + if self.reduction_config.subtract: + sq_q, sR_q, sdR_q, sFileName = self.loadRqz(self.reduction_config.subtract) + subtract = True + logging.warning(f'loaded background file: {sFileName}') + self.header.reduction.corrections.append(f'background from \'{sFileName}\' subtracted') + else: + subtract = False + + # load measurement data and do the reduction + datasetsRqz = [] + datasetsRlt = [] + for i, short_notation in enumerate(self.reduction_config.fileIdentifier): + logging.warning('reading input:') + self.header.measurement_data_files = [] + fromHDF.read_data(short_notation) + + if self.reduction_config.timeSlize: + wallTime_e = fromHDF.wallTime_e + columns = self.header.columns()+[fileio.Column('time', 's', 'time relative to start of measurement series')] + headerRqz = fileio.Orso(self.header.data_source(), self.header.reduction, columns) + + interval = self.reduction_config.timeSlize[0] + try: + start = self.reduction_config.timeSlize[1] + except: + start = 0 + try: + stop = self.reduction_config.timeSlize[2] + except: + stop = wallTime_e[-1] + # make overwriting log lines possible by removing newline at the end + logging.StreamHandler.terminator = "\r" + for i, time in enumerate(np.arange(start, stop, interval)): + logging.warning(f' time slize {i:4d}') + + filter_e = np.where((time0, filter_q, False) + if len(filter_q[filter_q]) > 0: + scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) + else: + logging.warning(f'# automatic scaling not possible') + scale = 1. + else: + filter_q = np.where(np.isnan(pR_q*R_q), False, True) + filter_q = np.where(R_q>0, filter_q, False) + filter_q = np.where(pR_q>0, filter_q, False) + if len(filter_q[filter_q]) > 0: + scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ + / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) + else: + logging.warning(f'# automatic scaling not possible') + scale = 1. + R_q /= scale + dR_q /= scale + logging.debug(f'# scaling factor = {scale}') + + return R_q, dR_q + + def project_on_qz(self, q_lz, R_lz, dR_lz, dq_lz, norm_lz, mask_lz): + q_q = self.grid.q() + mask_lzf = mask_lz.flatten() + q_lzf = q_lz.flatten()[mask_lzf] + R_lzf = R_lz.flatten()[mask_lzf] + dR_lzf = dR_lz.flatten()[mask_lzf] + dq_lzf = dq_lz.flatten()[mask_lzf] + norm_lzf = norm_lz.flatten()[mask_lzf] + + N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + + R_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf * R_lzf )[0] + R_q = R_q / N_q + + dR_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dR_lzf)**2 )[0] + dR_q = np.sqrt( dR_q ) / N_q + + # TODO: different error propagations for dR and dq! + N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf**2 )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + dq_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dq_lzf)**2 )[0] + dq_q = np.sqrt( dq_q / N_q ) + + q_q = 0.5 * (q_q + np.roll(q_q, 1)) + + return q_q[1:], R_q, dR_q, dq_q + + def loadRqz(self, name): + if os.path.exists(f'{self.data_reader_config.dataPath}/{name}'): + fileName = f'{self.data_reader_config.dataPath}/{name}' + elif os.path.exists(f'{self.data_reader_config.dataPath}/{name}.Rqz.ort'): + fileName = f'{self.data_reader_config.dataPath}/{name}.Rqz.ort' + else: + sys.exit(f'### the background file \'{self.data_reader_config.dataPath}/{name}\' does not exist! => stopping') + + q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) + + return q_q, Sq_q, dS_q, fileName + + def normalisation_map(self, short_notation): + dataPath = self.data_reader_config.dataPath + fromHDF = AmorData(self.startTime, self.header, self.data_reader_config) + normalisation_list = expand_file_list(short_notation) + name = str(normalisation_list[0]) + for i in range(1, len(normalisation_list), 1): + name = f'{name}_{normalisation_list[i]}' + if os.path.exists(f'{dataPath}/{name}.norm'): + logging.info(f'# normalisation matrix: found and using {dataPath}/{name}.norm') + norm_lz = np.loadtxt(f'{dataPath}/{name}.norm') + fh = open(f'{dataPath}/{name}.norm', 'r') + fh.readline() + normFileList = fh.readline().split('[')[1].split(']')[0].replace('\'', '').split(', ') + normAngle = float(fh.readline().split('= ')[1]) + fh.close() + for i, entry in enumerate(normFileList): + normFileList[i] = entry.split('/')[-1] + self.header.measurement_additional_files = normFileList + else: + logging.info(f'# normalisation matrix: using the files {normalisation_list}') + fromHDF.read_data(short_notation, norm=True) + normAngle = fromHDF.nu - fromHDF.mu + lamda_e = fromHDF.lamda_e + detZ_e = fromHDF.detZ_e + norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (self.grid.lamda(), self.grid.z())) + norm_lz = np.where(norm_lz>0, norm_lz, np.nan) + # correct for the SM reflectivity + lamda_l = self.grid.lamda() + theta_z = normAngle + fromHDF.delta_z + lamda_lz = (self.grid.lz().T*lamda_l[:-1]).T + theta_lz = self.grid.lz()*theta_z + qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz + 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) + norm_lz = norm_lz / Rsm_lz + if len(lamda_e) > 1e6: + head = ('normalisation matrix based on the measurements\n' + f'{fromHDF.file_list}\n' + f'nu - mu = {normAngle}\n' + f'shape= {np.shape(norm_lz)} (lambda, z)\n' + f'measured at mu = {fromHDF.mu:6.3f} deg\n' + f'N(l_lambda, z) = theta(z) / sum_i=-1..1 I(l_lambda+i, z)') + head = head.replace('../', '') + head = head.replace('./', '') + head = head.replace('raw/', '') + np.savetxt(f'{dataPath}/{name}.norm', norm_lz, header = head) + normFileList = fromHDF.file_list + return norm_lz, normAngle, normFileList + + + def project_on_lz(self, fromHDF, norm_lz, normAngle, lamda_e, detZ_e): + # projection on lambda-z-grid + lamda_l = self.grid.lamda() + theta_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z + lamda_lz = (self.grid.lz().T*lamda_l[:-1]).T + theta_lz = self.grid.lz()*theta_z + + thetaN_z = fromHDF.delta_z + normAngle + thetaN_lz = np.ones(np.shape(norm_lz))*thetaN_z + thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) + + mask_lz = np.where(np.isnan(norm_lz), False, True) + mask_lz = np.logical_and(mask_lz, np.where(np.absolute(thetaN_lz)>5e-3, True, False)) + mask_lz = np.logical_and(mask_lz, np.where(np.absolute(theta_lz)>5e-3, True, False)) + if self.reduction_config.thetaRange[1]<12: + mask_lz = np.logical_and(mask_lz, np.where(theta_lz >= self.reduction_config.thetaRange[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(theta_lz <= self.reduction_config.thetaRange[1], True, False)) + if self.reduction_config.thetaRangeR[1]<12: + t0 = fromHDF.nu - fromHDF.mu + mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 >= self.reduction_config.thetaRangeR[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 <= self.reduction_config.thetaRangeR[1], True, False)) + if self.reduction_config.lambdaRange[1]<15: + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= self.reduction_config.lambdaRange[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= self.reduction_config.lambdaRange[1], True, False)) + + # gravity correction + #theta_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) + theta_lz += np.rad2deg( np.arctan( 3.07e-10 * fromHDF.detectorDistance * lamda_lz**2 ) ) + + z_z = enumerate(theta_z) + qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz + int_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (lamda_l, self.grid.z())) + # cut normalisation sample horizon + int_lz = np.where(mask_lz, int_lz, np.nan) + thetaF_lz = np.where(mask_lz, theta_lz, np.nan) + + ref_lz = (int_lz * np.absolute(thetaN_lz)) / (norm_lz * np.absolute(thetaF_lz)) + err_lz = ref_lz * np.sqrt( 1/(int_lz+.1) + 1/norm_lz ) + + res_lz = np.ones((np.shape(lamda_l[:-1])[0], np.shape(theta_z)[0])) * 0.022**2 + res_lz = res_lz + (0.008/theta_lz)**2 + res_lz = qz_lz * np.sqrt(res_lz) + + return qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz - return R_q, dR_q diff --git a/neos.py b/neos.py index 85fb92b..57c5237 100644 --- a/neos.py +++ b/neos.py @@ -16,20 +16,9 @@ conventions (not strictly followed, yet): - to come """ - -import os -import sys -import logging -import logging.handlers -import numpy as np -from orsopy import fileio - -import libeos.reduction -from libeos.command_line import commandLineArgs, output_format_list -from libeos.dataset import AmorData, DataReaderConfig, Header -from libeos.instrument import Grid -from libeos.logconfig import setup_logging, update_loglevel -from libeos.reduction import autoscale, normalisation_map, project_on_lz, project_on_qz +from libeos.command_line import command_line_options +from libeos.logconfig import setup_logging +from libeos.reduction import AmorReduction #===================================================================================================== @@ -38,234 +27,16 @@ from libeos.reduction import autoscale, normalisation_map, project_on_lz, projec # - deal with background correction # - format of 'call' + add '-Y' if not supplied #===================================================================================================== -def loadRqz(name): - - if os.path.exists(f'{clas.dataPath}/{name}'): - fileName = f'{clas.dataPath}/{name}' - elif os.path.exists(f'{clas.dataPath}/{name}.Rqz.ort'): - fileName = f'{clas.dataPath}/{name}.Rqz.ort' - else: - sys.exit(f'### the background file \'{clas.dataPath}/{name}\' does not exist! => stopping') - - q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) - - return q_q, Sq_q, dS_q, fileName def main(): setup_logging() - global startTime, grid, clas, header - clas = commandLineArgs() - update_loglevel(clas.verbose, clas.debug) - - grid = Grid(clas.qResolution) - header = Header() - startTime = 0 - if not os.path.exists(f'{clas.dataPath}'): - os.system(f'mkdir {clas.dataPath}') - fromHDF = AmorData(startTime, header=header, config=DataReaderConfig( - year=clas.year, - dataPath=clas.dataPath, - sampleModel=clas.sampleModel, - chopperPhase=clas.chopperPhase, - chopperPhaseOffset=clas.chopperPhaseOffset, - yRange=clas.yRange, - lambdaRange=clas.lambdaRange, - qzRange=clas.qzRange, - offSpecular=clas.offSpecular, - mu=clas.mu, - nu=clas.nu, - muOffset=clas.muOffset - )) - logging.warning('\n######## eos - data reduction for Amor ########') - - # load or create normalisation matrix - if clas.normalisationFileIdentifier: - normalise = True - norm_lz, normAngle, normFileList = normalisation_map(clas.normalisationFileIdentifier[0], - header, grid, clas.dataPath) - header.reduction.corrections.append('normalisation with \'additional files\'') - else: - normalise = False - norm_lz = grid.lz() - normAngle = 1. - - logging.warning('normalisation matrix: none requested') - - # load R(q_z) curve to be subtracted: - if clas.subtract: - sq_q, sR_q, sdR_q, sFileName = loadRqz(clas.subtract) - subtract = True - logging.warning(f'loaded background file: {sFileName}') - header.reduction.corrections.append(f'background from \'{sFileName}\' subtracted') - else: - subtract = False - - # load measurement data and do the reduction - datasetsRqz = [] - datasetsRlt = [] - for i, short_notation in enumerate(clas.fileIdentifier): - logging.warning('reading input:') - header.measurement_data_files = [] - fromHDF.read_data(short_notation) - - if clas.timeSlize: - wallTime_e = fromHDF.wallTime_e - columns = header.columns() + [fileio.Column('time', 's', 'time relative to start of measurement series')] - headerRqz = fileio.Orso(header.data_source(), header.reduction, columns) - - interval = clas.timeSlize[0] - try: - start = clas.timeSlize[1] - except: - start = 0 - try: - stop = clas.timeSlize[2] - except: - stop = wallTime_e[-1] - # make overwriting log lines possible by removing newline at the end - logging.StreamHandler.terminator="\r" - for i, time in enumerate(np.arange(start, stop, interval)): - logging.warning(f' time slize {i:4d}') - - filter_e = np.where((time < wallTime_e) & (wallTime_e < time+interval), True, False) - lamda_e = fromHDF.lamda_e[filter_e] - detZ_e = fromHDF.detZ_e[filter_e] - - qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz = project_on_lz( - fromHDF, norm_lz, normAngle, lamda_e, detZ_e, - grid, clas.thetaRange, clas.thetaRangeR, clas.lambdaRange) - q_q, R_q, dR_q, dq_q = project_on_qz(qz_lz, ref_lz, err_lz, res_lz, norm_lz, mask_lz, grid) - - filter_q = np.where((clas.qzRange[0] < q_q) & (q_q < clas.qzRange[1]), True, False) - q_q = q_q[filter_q] - R_q = R_q[filter_q] - dR_q = dR_q[filter_q] - dq_q = dq_q[filter_q] - - if clas.autoscale: - R_q, dR_q = autoscale(q_q, R_q, dR_q, clas.autoscale) - - if subtract: - if len(q_q) == len(sq_q): - R_q -= sR_q - dR_q = np.sqrt( dR_q**2 + sdR_q**2 ) - else: - subtract = False - logging.warning(f'background file {sFileName} not compatible with q_z scale ({len(sq_q)} vs. {len(q_q)})') - - tme_q = np.ones(np.shape(q_q))*time - data = np.array([q_q, R_q, dR_q, dq_q, tme_q]).T - headerRqz.data_set = f'{i}: time = {time:8.1f} s to {time+interval:8.1f} s' - orso_data = fileio.OrsoDataset(headerRqz, data) - # make a copy of the header for the next iteration - headerRqz = fileio.Orso.from_dict(headerRqz.to_dict()) - datasetsRqz.append(orso_data) - # reset normal logging behavior - logging.StreamHandler.terminator="\n" - logging.warning(f' time slize {i:4d}, done') - else: - lamda_e = fromHDF.lamda_e - detZ_e = fromHDF.detZ_e - - qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz = project_on_lz( - fromHDF, norm_lz, normAngle, lamda_e, detZ_e, - grid, clas.thetaRange, clas.thetaRangeR, clas.lambdaRange - ) - try: - ref_lz *= clas.scale[i] - err_lz *= clas.scale[i] - except: - ref_lz *= clas.scale[-1] - err_lz *= clas.scale[-1] - - if 'Rqz.ort' in output_format_list(clas.outputFormat): - headerRqz = fileio.Orso(header.data_source(), header.reduction, header.columns()) - headerRqz.data_set = f'Nr {i} : mu = {fromHDF.mu:6.3f} deg' - - # projection on q-grid - q_q, R_q, dR_q, dq_q = project_on_qz(qz_lz, ref_lz, err_lz, res_lz, norm_lz, mask_lz, grid) - - filter_q = np.where((clas.qzRange[0] < q_q) & (q_q < clas.qzRange[1]), True, False) - q_q = q_q[filter_q] - R_q = R_q[filter_q] - dR_q = dR_q[filter_q] - dq_q = dq_q[filter_q] - - if libeos.reduction.autoscale: - if i == 0: - R_q, dR_q = autoscale(q_q, R_q, dR_q, clas.autoscale) - else: - pRq_z = datasetsRqz[i-1].data[:,1] - pdRq_z = datasetsRqz[i-1].data[:,2] - R_q, dR_q = autoscale(q_q, R_q, dR_q, clas.autoscale, pRq_z, pdRq_z) - - if subtract: - if len(q_q) == len(sq_q): - R_q -= sR_q - dR_q = np.sqrt( dR_q**2 + sdR_q**2 ) - else: - logging.warning(f'backgroung file {sFileName} not compatible with q_z scale ({len(sq_q)} vs. {len(q_q)})') - - data = np.array([q_q, R_q, dR_q, dq_q]).T - orso_data = fileio.OrsoDataset(headerRqz, data) - headerRqz = fileio.Orso(**headerRqz.to_dict()) - datasetsRqz.append(orso_data) - - if 'Rlt.ort' in output_format_list(clas.outputFormat): - columns = [ - fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), - fileio.Column('R', '', 'specular reflectivity'), - fileio.ErrorColumn(error_of='R', error_type='uncertainty', value_is='sigma'), - fileio.ErrorColumn(error_of='Qz', error_type='resolution', value_is='sigma'), - fileio.Column('lambda', 'angstrom', 'wavelength'), - fileio.Column('alpha_f', 'deg', 'final angle'), - fileio.Column('l', '', 'index of lambda-bin'), - fileio.Column('t', '', 'index of theta bin'), - fileio.Column('intensity', '', 'filtered neutron events per pixel'), - fileio.Column('norm', '', 'normalisation matrix'), - fileio.Column('mask', '', 'pixels used for calculating R(q_z)'), - ] - #data_source = fromHDF.data_source - headerRlt = fileio.Orso(header.data_source, header.reduction, columns) - - ts, zs=ref_lz.shape - lindex_lz=np.tile(np.arange(1, ts+1), (zs, 1)).T - tindex_lz=np.tile(np.arange(1, zs+1), (ts, 1)) - - j = 0 - for item in zip( - qz_lz.T, - ref_lz.T, - err_lz.T, - res_lz.T, - lamda_lz.T, - theta_lz.T, - lindex_lz.T, - tindex_lz.T, - int_lz.T, - norm_lz.T, - np.where(mask_lz, 1, 0).T - ): - data = np.array(list(item)).T - headerRlt = fileio.Orso(**headerRlt.to_dict()) - headerRlt.data_set = f'dataset_{i}_{j+1} : alpha_f = {theta_lz[0,j]:6.3f} deg' - orso_data = fileio.OrsoDataset(headerRlt, data) - datasetsRlt.append(orso_data) - j += 1 - - # output - logging.warning('output:') - - if 'Rqz.ort' in output_format_list(clas.outputFormat): - logging.warning(f' {clas.dataPath}/{clas.outputName}.Rqz.ort') - theSecondLine = f' {header.experiment.title} | {header.experiment.start_date} | sample {header.sample.name} | R(q_z)' - fileio.save_orso(datasetsRqz, f'{clas.dataPath}/{clas.outputName}.Rqz.ort', data_separator='\n', comment=theSecondLine) - - if 'Rlt.ort' in output_format_list(clas.outputFormat): - logging.warning(f' {clas.dataPath}/{clas.outputName}.Rlt.ort') - theSecondLine = f' {header.experiment.title} | {header.experiment.start_date} | sample {header.sample.name} | R(lambda, theta)' - fileio.save_orso(datasetsRlt, f'{clas.dataPath}/{clas.outputName}.Rlt.ort', data_separator='\n', comment=theSecondLine) + # read command line arguments and generate classes holding configuration parameters + reader_config, reduction_config, output_config = command_line_options() + # Create reducer with these arguments + reducer = AmorReduction(reader_config, reduction_config, output_config) + # Perform actual reduction + reducer.reduce() if __name__ == '__main__': main() From a56b10ae9a159f44726589f50a5580cf4c9ce9e8 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 13:48:13 +0100 Subject: [PATCH 06/17] Separate reader from experiment config --- libeos/command_line.py | 9 +++---- libeos/dataset.py | 13 +++++----- libeos/options.py | 17 +++++++++---- libeos/reduction.py | 54 ++++++++++++++++++++---------------------- neos.py | 8 +++++-- 5 files changed, 57 insertions(+), 44 deletions(-) diff --git a/libeos/command_line.py b/libeos/command_line.py index d43053d..de8bc92 100644 --- a/libeos/command_line.py +++ b/libeos/command_line.py @@ -2,7 +2,7 @@ import argparse from datetime import datetime from .logconfig import update_loglevel -from .options import DataReaderConfig, OutputConfig, ReductionConfig +from .options import ReaderConfig, EOSConfig, ExperimentConfig, OutputConfig, ReductionConfig def commandLineArgs(): @@ -160,9 +160,10 @@ def command_line_options(): clas = commandLineArgs() update_loglevel(clas.verbose, clas.debug) - reader_config = DataReaderConfig( + reader_config = ReaderConfig( year=clas.year, - dataPath=clas.dataPath, + dataPath=clas.dataPath) + experiment_config = ExperimentConfig( sampleModel=clas.sampleModel, chopperPhase=clas.chopperPhase, chopperPhaseOffset=clas.chopperPhaseOffset, @@ -190,4 +191,4 @@ def command_line_options(): outputFormats=output_format_list(clas.outputFormat), outputName=clas.outputName ) - return reader_config, reduction_config, output_config \ No newline at end of file + return EOSConfig(reader_config, experiment_config, reduction_config, output_config) \ No newline at end of file diff --git a/libeos/dataset.py b/libeos/dataset.py index 10cf34d..6c06fe4 100644 --- a/libeos/dataset.py +++ b/libeos/dataset.py @@ -11,7 +11,7 @@ from orsopy import fileio from . import __version__, const from .instrument import Detector -from .options import DataReaderConfig +from .options import ExperimentConfig, ReaderConfig class Header: @@ -62,10 +62,11 @@ class Header: class AmorData: """read meta-data and event streams from .hdf file(s), apply filters and conversions""" #------------------------------------------------------------------------------------------------- - def __init__(self, startTime, header:Header, config: DataReaderConfig): + def __init__(self, startTime, header:Header, reader_config: ReaderConfig, config: ExperimentConfig): self.startTime = startTime self.header = header self.config = config + self.reader_config = reader_config #------------------------------------------------------------------------------------------------- def read_data(self, short_notation, norm=False): self.data_file_numbers = self.expand_file_list(short_notation) @@ -93,16 +94,16 @@ class AmorData: #------------------------------------------------------------------------------------------------- def path_generator(self, number): - fileName = f'amor{self.config.year}n{number:06d}.hdf' - if os.path.exists(f'{self.config.dataPath}/{fileName}'): - path = self.config.dataPath + fileName = f'amor{self.reader_config.year}n{number:06d}.hdf' + if os.path.exists(f'{self.reader_config.dataPath}/{fileName}'): + path = self.reader_config.dataPath elif os.path.exists(fileName): path = '.' elif os.path.exists(f'./raw/{fileName}'): path = './raw' elif os.path.exists(f'../raw/{fileName}'): path = '../raw' - elif os.path.exists(f'/afs/psi.ch/project/sinqdata/{self.config.year}/amor/{int(number/1000)}/{fileName}'): + elif os.path.exists(f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}/{fileName}'): path = '/afs/psi.ch/project/sinqdata/{self.config.year}/amor/{int(number/1000)}' else: sys.exit(f'# ERROR: the file {fileName} is nowhere to be found!') diff --git a/libeos/options.py b/libeos/options.py index 606d6da..518ab1d 100644 --- a/libeos/options.py +++ b/libeos/options.py @@ -4,11 +4,13 @@ Classes for stroing various configurations needed for reduction. from dataclasses import dataclass, field from typing import Optional, Tuple - @dataclass -class DataReaderConfig: +class ReaderConfig: year: int dataPath: str + +@dataclass +class ExperimentConfig: sampleModel: str chopperPhase: float @@ -31,7 +33,7 @@ class ReductionConfig: lambdaRange: Tuple[float, float] fileIdentifier: list = field(default_factory=lambda: ["0"]) - scale: list = field(default_factory=lambda: [1]) #per file scaling, if less elements then files use the last one + scale: list = field(default_factory=lambda: [1]) #per file scaling; if less elements than files use the last one subtract: Optional[str] = None normalisationFileIdentifier: Optional[list] = None @@ -39,5 +41,12 @@ class ReductionConfig: @dataclass class OutputConfig: - outputFormats: list #output_format_list(clas.outputFormat) + outputFormats: list outputName: str + +@dataclass +class EOSConfig: + reader: ReaderConfig + experiment: ExperimentConfig + reductoin: ReductionConfig + output: OutputConfig \ No newline at end of file diff --git a/libeos/reduction.py b/libeos/reduction.py index e41a526..1f62511 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -7,26 +7,24 @@ from orsopy import fileio from .command_line import expand_file_list from .dataset import AmorData, Header -from .options import DataReaderConfig, ReductionConfig, OutputConfig +from .options import EOSConfig from .instrument import Grid class AmorReduction: - def __init__(self, - data_reader_config: DataReaderConfig, - reduction_config: ReductionConfig, - output_config: OutputConfig): - self.data_reader_config = data_reader_config - self.reduction_config = reduction_config - self.output_config = output_config - - def reduce(self): - self.grid = Grid(self.reduction_config.qResolution) + def __init__(self, config: EOSConfig): + self.experiment_config = config.experiment + self.reader_config = config.reader + self.reduction_config = config.reductoin + self.output_config = config.output + self.grid = Grid(config.reductoin.qResolution) self.header = Header() self.startTime = 0 - if not os.path.exists(f'{self.data_reader_config.dataPath}'): - os.system(f'mkdir {self.data_reader_config.dataPath}') - fromHDF = AmorData(self.startTime, header=self.header, config=self.data_reader_config) - logging.warning('\n######## eos - data reduction for Amor ########') + + def reduce(self): + if not os.path.exists(f'{self.reader_config.dataPath}'): + logging.debug(f'Creating destination path {self.reader_config.dataPath}') + os.system(f'mkdir {self.reader_config.dataPath}') + fromHDF = AmorData(self.startTime, header=self.header, reader_config=self.reader_config, config=self.experiment_config) # load or create normalisation matrix if self.reduction_config.normalisationFileIdentifier: @@ -84,7 +82,7 @@ class AmorReduction: fromHDF, norm_lz, normAngle, lamda_e, detZ_e) q_q, R_q, dR_q, dq_q = self.project_on_qz(qz_lz, ref_lz, err_lz, res_lz, norm_lz, mask_lz) - filter_q = np.where((self.data_reader_config.qzRange[0] stopping') + sys.exit(f'### the background file \'{self.reader_config.dataPath}/{name}\' does not exist! => stopping') q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) return q_q, Sq_q, dS_q, fileName def normalisation_map(self, short_notation): - dataPath = self.data_reader_config.dataPath - fromHDF = AmorData(self.startTime, self.header, self.data_reader_config) + dataPath = self.reader_config.dataPath + fromHDF = AmorData(self.startTime, header=self.header, reader_config=self.reader_config, config=self.experiment_config) normalisation_list = expand_file_list(short_notation) name = str(normalisation_list[0]) for i in range(1, len(normalisation_list), 1): diff --git a/neos.py b/neos.py index 57c5237..8264173 100644 --- a/neos.py +++ b/neos.py @@ -15,6 +15,7 @@ conventions (not strictly followed, yet): _q = q_z - to come """ +import logging from libeos.command_line import command_line_options from libeos.logconfig import setup_logging @@ -30,13 +31,16 @@ from libeos.reduction import AmorReduction def main(): setup_logging() + logging.warning('######## eos - data reduction for Amor ########') # read command line arguments and generate classes holding configuration parameters - reader_config, reduction_config, output_config = command_line_options() + config = command_line_options() # Create reducer with these arguments - reducer = AmorReduction(reader_config, reduction_config, output_config) + reducer = AmorReduction(config) # Perform actual reduction reducer.reduce() + logging.info('######## eos - finished ########') + if __name__ == '__main__': main() From 44ab8fb211534454e02389b57c12fbc1405c7c3f Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 15:26:29 +0100 Subject: [PATCH 07/17] Start refactoring reduction into multiple methods --- libeos/dataset.py | 10 +- libeos/reduction.py | 380 ++++++++++++++++++++++---------------------- 2 files changed, 201 insertions(+), 189 deletions(-) diff --git a/libeos/dataset.py b/libeos/dataset.py index 6c06fe4..117f509 100644 --- a/libeos/dataset.py +++ b/libeos/dataset.py @@ -13,7 +13,6 @@ from . import __version__, const from .instrument import Detector from .options import ExperimentConfig, ReaderConfig - class Header: """orso compatible output file header content""" @@ -58,6 +57,15 @@ class Header: ] return cols + def orso_header(self, columns=None, extra_columns=[]): + """ + Generate ORSO header from a copy of this class' data. + """ + ds = fileio.DataSource.from_dict(self.data_source().to_dict()) + red = fileio.Reduction.from_dict(self.reduction.to_dict()) + if columns is None: + columns = self.columns() + return fileio.Orso(ds, red, columns+extra_columns) class AmorData: """read meta-data and event streams from .hdf file(s), apply filters and conversions""" diff --git a/libeos/reduction.py b/libeos/reduction.py index 1f62511..938e78c 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -24,197 +24,202 @@ class AmorReduction: if not os.path.exists(f'{self.reader_config.dataPath}'): logging.debug(f'Creating destination path {self.reader_config.dataPath}') os.system(f'mkdir {self.reader_config.dataPath}') - fromHDF = AmorData(self.startTime, header=self.header, reader_config=self.reader_config, config=self.experiment_config) # load or create normalisation matrix if self.reduction_config.normalisationFileIdentifier: - normalise = True - norm_lz, normAngle, normFileList = self.normalisation_map(self.reduction_config.normalisationFileIdentifier[0]) - self.header.reduction.corrections.append('normalisation with \'additional files\'') + self.create_normalisation_map(self.reduction_config.normalisationFileIdentifier[0]) else: - normalise = False - norm_lz = self.grid.lz() - normAngle = 1. + self.norm_lz = self.grid.lz() + self.normAngle = 1. logging.warning('normalisation matrix: none requested') # load R(q_z) curve to be subtracted: if self.reduction_config.subtract: - sq_q, sR_q, sdR_q, sFileName = self.loadRqz(self.reduction_config.subtract) - subtract = True - logging.warning(f'loaded background file: {sFileName}') - self.header.reduction.corrections.append(f'background from \'{sFileName}\' subtracted') + self.sq_q, self.sR_q, self.sdR_q, self.sFileName = self.loadRqz(self.reduction_config.subtract) + logging.warning(f'loaded background file: {self.sFileName}') + self.header.reduction.corrections.append(f'background from \'{self.sFileName}\' subtracted') + self.subtract = True else: - subtract = False + self.subtract = False # load measurement data and do the reduction - datasetsRqz = [] - datasetsRlt = [] + self.datasetsRqz = [] + self.datasetsRlt = [] + self.file_reader = AmorData(self.startTime, header=self.header, reader_config=self.reader_config, config=self.experiment_config) for i, short_notation in enumerate(self.reduction_config.fileIdentifier): - logging.warning('reading input:') - self.header.measurement_data_files = [] - fromHDF.read_data(short_notation) - - if self.reduction_config.timeSlize: - wallTime_e = fromHDF.wallTime_e - columns = self.header.columns()+[fileio.Column('time', 's', 'time relative to start of measurement series')] - headerRqz = fileio.Orso(self.header.data_source(), self.header.reduction, columns) - - interval = self.reduction_config.timeSlize[0] - try: - start = self.reduction_config.timeSlize[1] - except: - start = 0 - try: - stop = self.reduction_config.timeSlize[2] - except: - stop = wallTime_e[-1] - # make overwriting log lines possible by removing newline at the end - logging.StreamHandler.terminator = "\r" - for i, time in enumerate(np.arange(start, stop, interval)): - logging.warning(f' time slize {i:4d}') - - filter_e = np.where((time0, norm_lz, np.nan) + self.norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (self.grid.lamda(), self.grid.z())) + self.norm_lz = np.where(self.norm_lz>0, self.norm_lz, np.nan) # correct for the SM reflectivity lamda_l = self.grid.lamda() - theta_z = normAngle + fromHDF.delta_z + theta_z = self.normAngle + fromHDF.delta_z lamda_lz = (self.grid.lz().T*lamda_l[:-1]).T theta_lz = self.grid.lz()*theta_z qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz 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) - norm_lz = norm_lz / Rsm_lz + self.norm_lz = self.norm_lz / Rsm_lz if len(lamda_e) > 1e6: head = ('normalisation matrix based on the measurements\n' f'{fromHDF.file_list}\n' - f'nu - mu = {normAngle}\n' - f'shape= {np.shape(norm_lz)} (lambda, z)\n' + f'nu - mu = {self.normAngle}\n' + f'shape= {np.shape(self.norm_lz)} (lambda, z)\n' f'measured at mu = {fromHDF.mu:6.3f} deg\n' f'N(l_lambda, z) = theta(z) / sum_i=-1..1 I(l_lambda+i, z)') head = head.replace('../', '') head = head.replace('./', '') head = head.replace('raw/', '') - np.savetxt(f'{dataPath}/{name}.norm', norm_lz, header = head) - normFileList = fromHDF.file_list - return norm_lz, normAngle, normFileList - + np.savetxt(f'{dataPath}/{name}.norm', self.norm_lz, header = head) + self.normFileList = fromHDF.file_list + self.header.reduction.corrections.append('normalisation with \'additional files\'') def project_on_lz(self, fromHDF, norm_lz, normAngle, lamda_e, detZ_e): # projection on lambda-z-grid From 6cee00adaa30b9674f95886a5312d30e054978c6 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 16:39:30 +0100 Subject: [PATCH 08/17] Move startTime to ReaderConfig, separte Header into module, AmorData class to read directly (clearify scope) --- libeos/dataset.py | 97 ++++++++++---------------------------------- libeos/header.py | 66 ++++++++++++++++++++++++++++++ libeos/instrument.py | 11 +++-- libeos/options.py | 1 + libeos/reduction.py | 20 +++++---- 5 files changed, 106 insertions(+), 89 deletions(-) create mode 100644 libeos/header.py diff --git a/libeos/dataset.py b/libeos/dataset.py index 117f509..fd15377 100644 --- a/libeos/dataset.py +++ b/libeos/dataset.py @@ -1,6 +1,5 @@ import logging import os -import platform import subprocess import sys from datetime import datetime @@ -9,76 +8,25 @@ import h5py import numpy as np from orsopy import fileio -from . import __version__, const +from . import const +from .header import Header from .instrument import Detector from .options import ExperimentConfig, ReaderConfig -class Header: - """orso compatible output file header content""" - - def __init__(self): - self.owner = None - self.experiment = None - self.sample = None - self.measurement_instrument_settings = None - self.measurement_scheme = None - self.measurement_data_files = [] - self.measurement_additional_files = [] - - self.reduction = fileio.Reduction( - software = fileio.Software('eos', version=__version__), - call = ' '.join(sys.argv), - computer = platform.node(), - timestamp = datetime.now(), - creator = None, - corrections = ['histogramming in lambda and alpha_f', - 'gravity'], - ) - #------------------------------------------------------------------------------------------------- - def data_source(self): - return fileio.DataSource( - self.owner, - self.experiment, - self.sample, - fileio.Measurement( - instrument_settings = self.measurement_instrument_settings, - scheme = self.measurement_scheme, - data_files = self.measurement_data_files, - additional_files = self.measurement_additional_files, - ), - ) - #------------------------------------------------------------------------------------------------- - def columns(self): - cols = [ - fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), - fileio.Column('R', '', 'specular reflectivity'), - fileio.ErrorColumn(error_of='R', error_type='uncertainty', distribution='gaussian', value_is='sigma'), - fileio.ErrorColumn(error_of='Qz', error_type='resolution', distribution='gaussian', value_is='sigma'), - ] - return cols - - def orso_header(self, columns=None, extra_columns=[]): - """ - Generate ORSO header from a copy of this class' data. - """ - ds = fileio.DataSource.from_dict(self.data_source().to_dict()) - red = fileio.Reduction.from_dict(self.reduction.to_dict()) - if columns is None: - columns = self.columns() - return fileio.Orso(ds, red, columns+extra_columns) class AmorData: """read meta-data and event streams from .hdf file(s), apply filters and conversions""" #------------------------------------------------------------------------------------------------- - def __init__(self, startTime, header:Header, reader_config: ReaderConfig, config: ExperimentConfig): - self.startTime = startTime + def __init__(self, header: Header, reader_config: ReaderConfig, config: ExperimentConfig, short_notation, norm=False): + self.startTime = reader_config.startTime self.header = header self.config = config self.reader_config = reader_config + self.read_data(short_notation, norm=norm) + #------------------------------------------------------------------------------------------------- def read_data(self, short_notation, norm=False): self.data_file_numbers = self.expand_file_list(short_notation) - #self.year = clas.year self.file_list = [] for number in self.data_file_numbers: self.file_list.append(self.path_generator(number)) @@ -107,15 +55,15 @@ class AmorData: path = self.reader_config.dataPath elif os.path.exists(fileName): path = '.' - elif os.path.exists(f'./raw/{fileName}'): - path = './raw' - elif os.path.exists(f'../raw/{fileName}'): - path = '../raw' + elif os.path.exists(os.path.join('.','raw', fileName)): + path = os.path.join('.','raw') + elif os.path.exists(os.path.join('..','raw', fileName)): + path = os.path.join('..','raw') elif os.path.exists(f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}/{fileName}'): - path = '/afs/psi.ch/project/sinqdata/{self.config.year}/amor/{int(number/1000)}' + path = f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}' else: sys.exit(f'# ERROR: the file {fileName} is nowhere to be found!') - return f'{path}/{fileName}' + return os.path.join(path, fileName) #------------------------------------------------------------------------------------------------- def expand_file_list(self, short_notation): """Evaluate string entry for file number lists""" @@ -135,17 +83,16 @@ class AmorData: #------------------------------------------------------------------------------------------------- def resolve_pixels(self): """determine spatial coordinats and angles from pixel number""" - det = Detector() - nPixel = det.nWires * det.nStripes * det.nBlades + nPixel = Detector.nWires * Detector.nStripes * Detector.nBlades pixelID = np.arange(nPixel) - (bladeNr, bPixel) = np.divmod(pixelID, det.nWires * det.nStripes) - (bZi, detYi) = np.divmod(bPixel, det.nStripes) # z index on blade, y index on detector - detZi = bladeNr * det.nWires + bZi # z index on detector - detX = bZi * det.dX # x position in detector - # detZ = det.zero - bladeNr * det.bladeZ - bZi * det.dZ # z position on detector - bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*det.bladeZ / det.distance) ) - delta = (det.nBlades/2. - bladeNr) * bladeAngle \ - - np.rad2deg( np.arctan(bZi*det.dZ / ( det.distance + bZi * det.dX) ) ) + (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] return np.vstack((detYi.T, detZi.T, detX.T, delta.T)).T #return matr @@ -268,7 +215,7 @@ class AmorData: # TODO: should extract monitor from counts or beam current times time self.monitor1 = self.ctime - self.monitor2=self.monitor1 + self.monitor2 = self.monitor1 # read data event streams tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 diff --git a/libeos/header.py b/libeos/header.py new file mode 100644 index 0000000..527f070 --- /dev/null +++ b/libeos/header.py @@ -0,0 +1,66 @@ +""" +Class to handle Orso header information that changes gradually during the reduction process. +""" + +import platform +import sys +from datetime import datetime + +from orsopy import fileio + +from libeos import __version__ + + +class Header: + """orso compatible output file header content""" + + def __init__(self): + self.owner = None + self.experiment = None + self.sample = None + self.measurement_instrument_settings = None + self.measurement_scheme = None + self.measurement_data_files = [] + self.measurement_additional_files = [] + + self.reduction = fileio.Reduction( + software = fileio.Software('eos', version=__version__), + call = ' '.join(sys.argv), + computer = platform.node(), + timestamp = datetime.now(), + creator = None, + corrections = ['histogramming in lambda and alpha_f', + 'gravity'], + ) + #------------------------------------------------------------------------------------------------- + def data_source(self): + return fileio.DataSource( + self.owner, + self.experiment, + self.sample, + fileio.Measurement( + instrument_settings = self.measurement_instrument_settings, + scheme = self.measurement_scheme, + data_files = self.measurement_data_files, + additional_files = self.measurement_additional_files, + ), + ) + #------------------------------------------------------------------------------------------------- + def columns(self): + cols = [ + fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), + fileio.Column('R', '', 'specular reflectivity'), + fileio.ErrorColumn(error_of='R', error_type='uncertainty', distribution='gaussian', value_is='sigma'), + fileio.ErrorColumn(error_of='Qz', error_type='resolution', distribution='gaussian', value_is='sigma'), + ] + return cols + + def orso_header(self, columns=None, extra_columns=[]): + """ + Generate ORSO header from a copy of this class' data. + """ + ds = fileio.DataSource.from_dict(self.data_source().to_dict()) + red = fileio.Reduction.from_dict(self.reduction.to_dict()) + if columns is None: + columns = self.columns() + return fileio.Orso(ds, red, columns+extra_columns) diff --git a/libeos/instrument.py b/libeos/instrument.py index 7bd2c77..9a48817 100644 --- a/libeos/instrument.py +++ b/libeos/instrument.py @@ -21,7 +21,6 @@ class Detector: class Grid: def __init__(self, qResolution): - self.det = Detector() self.lamdaCut = const.lamdaCut self.dldl = 0.005 # Delta lambda / lambda self.qResolution = qResolution @@ -46,24 +45,24 @@ class Grid: return lamda_grid def z(self): - return np.arange(self.det.nBlades*self.det.nWires+1) + return np.arange(Detector.nBlades*Detector.nWires+1) def lz(self): return np.ones(( np.shape(self.lamda()[:-1])[0], np.shape(self.z()[:-1])[0] )) def delta(self, detectorDistance): # unused for now - bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*self.det.bladeZ / detectorDistance) ) - blade_grid = np.arctan( np.arange(33) * self.det.dZ / ( detectorDistance + np.arange(33) * self.det.dX) ) + bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*Detector.bladeZ / detectorDistance) ) + blade_grid = np.arctan( np.arange(33) * Detector.dZ / ( detectorDistance + np.arange(33) * Detector.dX) ) blade_grid = np.rad2deg(blade_grid) stepWidth = blade_grid[1] - blade_grid[0] blade_grid = blade_grid - 0.2 * stepWidth delta_grid = [] - for b in np.arange(self.det.nBlades-1): + for b in np.arange(Detector.nBlades-1): delta_grid = np.concatenate((delta_grid, blade_grid), axis=None) blade_grid = blade_grid + bladeAngle delta_grid = delta_grid[delta_grid Date: Mon, 4 Mar 2024 16:55:51 +0100 Subject: [PATCH 09/17] Rename to file_reader and add annotations to AmorData --- libeos/{dataset.py => file_reader.py} | 39 ++++++++++++++++++++++----- libeos/header.py | 2 +- libeos/reduction.py | 2 +- 3 files changed, 34 insertions(+), 9 deletions(-) rename libeos/{dataset.py => file_reader.py} (94%) diff --git a/libeos/dataset.py b/libeos/file_reader.py similarity index 94% rename from libeos/dataset.py rename to libeos/file_reader.py index fd15377..edc7922 100644 --- a/libeos/dataset.py +++ b/libeos/file_reader.py @@ -3,6 +3,7 @@ import os import subprocess import sys from datetime import datetime +from typing import List import h5py import numpy as np @@ -16,17 +17,41 @@ from .options import ExperimentConfig, ReaderConfig class AmorData: """read meta-data and event streams from .hdf file(s), apply filters and conversions""" + chopperDetectorDistance: float + chopperDistance: float + chopperPhase: float + chopperSpeed: float + ctime: 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 + monitor1: float + monitor2: float + mu: float + nu: float + tau: float + tofCut: float + start_date: str + #------------------------------------------------------------------------------------------------- - def __init__(self, header: Header, reader_config: ReaderConfig, config: ExperimentConfig, short_notation, norm=False): + 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 - self.read_data(short_notation, norm=norm) + self.expand_file_list(short_notation) + self.read_data(norm=norm) #------------------------------------------------------------------------------------------------- - def read_data(self, short_notation, norm=False): - self.data_file_numbers = self.expand_file_list(short_notation) + def read_data(self, norm=False): self.file_list = [] for number in self.data_file_numbers: self.file_list.append(self.path_generator(number)) @@ -79,7 +104,7 @@ class AmorData: file_list += range(int(i.split('-', 1)[0]), int(i.split('-', 1)[1])+1, int(step)) else: file_list += [int(i)] - return sorted(file_list) + self.data_file_numbers=sorted(file_list) #------------------------------------------------------------------------------------------------- def resolve_pixels(self): """determine spatial coordinats and angles from pixel number""" @@ -120,7 +145,7 @@ class AmorData: sourceProbe = 'neutron' start_time = self.hdf['entry1/start_time'][0].decode('utf-8') - start_date = start_time.split(' ')[0] + self.start_date = start_time.split(' ')[0] if self.config.sampleModel: model = self.config.sampleModel @@ -138,7 +163,7 @@ class AmorData: self.header.experiment = fileio.Experiment( title = title, instrument = instrumentName, - start_date = start_date, + start_date = self.start_date, probe = sourceProbe, facility = source, proposalID = proposal_id diff --git a/libeos/header.py b/libeos/header.py index 527f070..e46fc5e 100644 --- a/libeos/header.py +++ b/libeos/header.py @@ -8,7 +8,7 @@ from datetime import datetime from orsopy import fileio -from libeos import __version__ +from . import __version__ class Header: diff --git a/libeos/reduction.py b/libeos/reduction.py index 54245e6..321ccee 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -6,7 +6,7 @@ import numpy as np from orsopy import fileio from .command_line import expand_file_list -from .dataset import AmorData +from .file_reader import AmorData from .header import Header from .options import EOSConfig from .instrument import Grid From 8300e842dc32adf828284715f1aaeb2a70c74551 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 17:39:10 +0100 Subject: [PATCH 10/17] Add test cases for full reduction with and w/o slicing (data to be stored in "test_data" folder) --- .gitignore | 1 + libeos/command_line.py | 1 - libeos/file_reader.py | 4 +- libeos/options.py | 6 +-- libeos/reduction.py | 45 ++++++++++---------- tests/__init__.py | 0 tests/test_full_analysis.py | 82 +++++++++++++++++++++++++++++++++++++ 7 files changed, 108 insertions(+), 31 deletions(-) create mode 100644 tests/__init__.py create mode 100644 tests/test_full_analysis.py diff --git a/.gitignore b/.gitignore index 0fbb9a8..cfe693c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ *~ *.bak +*.log __py_cache__ raw .idea diff --git a/libeos/command_line.py b/libeos/command_line.py index de8bc92..7fae9e5 100644 --- a/libeos/command_line.py +++ b/libeos/command_line.py @@ -180,7 +180,6 @@ def command_line_options(): autoscale=clas.autoscale, thetaRange=clas.thetaRange, thetaRangeR=clas.thetaRangeR, - lambdaRange=clas.lambdaRange, fileIdentifier=clas.fileIdentifier, scale=clas.scale, subtract=clas.subtract, diff --git a/libeos/file_reader.py b/libeos/file_reader.py index edc7922..e9f49ca 100644 --- a/libeos/file_reader.py +++ b/libeos/file_reader.py @@ -76,7 +76,7 @@ class AmorData: #------------------------------------------------------------------------------------------------- def path_generator(self, number): fileName = f'amor{self.reader_config.year}n{number:06d}.hdf' - if os.path.exists(f'{self.reader_config.dataPath}/{fileName}'): + if os.path.exists(os.path.join(self.reader_config.dataPath,fileName)): path = self.reader_config.dataPath elif os.path.exists(fileName): path = '.' @@ -149,8 +149,6 @@ class AmorData: if self.config.sampleModel: model = self.config.sampleModel - else: - model = None # assembling orso header information self.header.owner = fileio.Person( diff --git a/libeos/options.py b/libeos/options.py index 022eee0..22c9208 100644 --- a/libeos/options.py +++ b/libeos/options.py @@ -12,13 +12,12 @@ class ReaderConfig: @dataclass class ExperimentConfig: - sampleModel: str - chopperPhase: float yRange: Tuple[float, float] lambdaRange: Tuple[float, float] qzRange: Tuple[float, float] + sampleModel: Optional[str] = None chopperPhaseOffset: float = 0.0 mu: Optional[float] = None nu: Optional[float] = None @@ -28,14 +27,13 @@ class ExperimentConfig: @dataclass class ReductionConfig: qResolution: float - autoscale: Tuple[bool, bool] thetaRange: Tuple[float, float] thetaRangeR: Tuple[float, float] - lambdaRange: Tuple[float, float] fileIdentifier: list = field(default_factory=lambda: ["0"]) scale: list = field(default_factory=lambda: [1]) #per file scaling; if less elements than files use the last one + autoscale: Optional[Tuple[bool, bool]] = None subtract: Optional[str] = None normalisationFileIdentifier: Optional[list] = None timeSlize: Optional[list] = None diff --git a/libeos/reduction.py b/libeos/reduction.py index 321ccee..44c6073 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -210,18 +210,16 @@ class AmorReduction: logging.warning(f' time slize {ti:4d}, done') def save_Rqz(self): - logging.warning(f' {self.reader_config.dataPath}/{self.output_config.outputName}.Rqz.ort') + fname = os.path.join(self.reader_config.dataPath, f'{self.output_config.outputName}.Rqz.ort') + logging.warning(f' {fname}') theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | R(q_z)' - fileio.save_orso(self.datasetsRqz, f'{self.reader_config.dataPath}/{self.output_config.outputName}.Rqz.ort', - data_separator='\n', - comment=theSecondLine) + fileio.save_orso(self.datasetsRqz, fname, data_separator='\n', comment=theSecondLine) def save_Rtl(self): - logging.warning(f' {self.reader_config.dataPath}/{self.output_config.outputName}.Rlt.ort') + fname = os.path.join(self.reader_config.dataPath, f'{self.output_config.outputName}.Rlt.ort') + logging.warning(f' {fname}') theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | R(lambda, theta)' - fileio.save_orso(self.datasetsRlt, f'{self.reader_config.dataPath}/{self.output_config.outputName}.Rlt.ort', - data_separator='\n', - comment=theSecondLine) + fileio.save_orso(self.datasetsRlt, fname, data_separator='\n', comment=theSecondLine) def autoscale(self, q_q, R_q, dR_q, pR_q=[], pdR_q=[]): autoscale = self.reduction_config.autoscale @@ -278,12 +276,13 @@ class AmorReduction: return q_q[1:], R_q, dR_q, dq_q def loadRqz(self, name): - if os.path.exists(f'{self.reader_config.dataPath}/{name}'): - fileName = f'{self.reader_config.dataPath}/{name}' - elif os.path.exists(f'{self.reader_config.dataPath}/{name}.Rqz.ort'): - fileName = f'{self.reader_config.dataPath}/{name}.Rqz.ort' + fname = os.path.join(self.reader_config.dataPath, name) + if os.path.exists(fname): + fileName = fname + elif os.path.exists(f'{fname}.Rqz.ort'): + fileName = f'{fname}.Rqz.ort' else: - sys.exit(f'### the background file \'{self.reader_config.dataPath}/{name}\' does not exist! => stopping') + sys.exit(f'### the background file \'{fname}\' does not exist! => stopping') q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) @@ -295,14 +294,14 @@ class AmorReduction: name = str(normalisation_list[0]) for i in range(1, len(normalisation_list), 1): name = f'{name}_{normalisation_list[i]}' - if os.path.exists(f'{dataPath}/{name}.norm'): - logging.info(f'# normalisation matrix: found and using {dataPath}/{name}.norm') + n_path = os.path.join(dataPath, f'{name}.norm') + if os.path.exists(n_path): + logging.info(f'# normalisation matrix: found and using {n_path}') self.norm_lz = np.loadtxt(f'{dataPath}/{name}.norm') - fh = open(f'{dataPath}/{name}.norm', 'r') - fh.readline() - self.normFileList = fh.readline().split('[')[1].split(']')[0].replace('\'', '').split(', ') - self.normAngle = float(fh.readline().split('= ')[1]) - fh.close() + with open(n_path, 'r') as fh: + fh.readline() + self.normFileList = fh.readline().split('[')[1].split(']')[0].replace('\'', '').split(', ') + self.normAngle = float(fh.readline().split('= ')[1]) for i, entry in enumerate(self.normFileList): self.normFileList[i] = entry.split('/')[-1] self.header.measurement_additional_files = self.normFileList @@ -362,9 +361,9 @@ class AmorReduction: t0 = fromHDF.nu - fromHDF.mu mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 >= self.reduction_config.thetaRangeR[0], True, False)) mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 <= self.reduction_config.thetaRangeR[1], True, False)) - if self.reduction_config.lambdaRange[1]<15: - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= self.reduction_config.lambdaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= self.reduction_config.lambdaRange[1], True, False)) + if self.experiment_config.lambdaRange[1]<15: + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= self.experiment_config.lambdaRange[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= self.experiment_config.lambdaRange[1], True, False)) # gravity correction #theta_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_full_analysis.py b/tests/test_full_analysis.py new file mode 100644 index 0000000..29606c9 --- /dev/null +++ b/tests/test_full_analysis.py @@ -0,0 +1,82 @@ +import os +from unittest import TestCase +from libeos import options, reduction, logconfig + +logconfig.setup_logging() +logconfig.update_loglevel(True, False) + +class FullAmorTest(TestCase): + def setUp(self): + self.reader_config = options.ReaderConfig( + year=2023, + dataPath=os.path.join('..', "test_data")) + + def tearDown(self): + for fi in ['test.Rqz.ort', '614.norm']: + try: + os.unlink(os.path.join(self.reader_config.dataPath, fi)) + except FileNotFoundError: + pass + + + def test_time_slicing(self): + experiment_config = options.ExperimentConfig( + chopperPhase=-13.5, + chopperPhaseOffset=-5, + yRange=(11., 41.), + lambdaRange=(2., 15.), + qzRange=(0.005, 0.30), + offSpecular=False, + mu=0, + nu=0, + muOffset=0.0, + sampleModel='air | 10 H2O | D2O' + ) + reduction_config = options.ReductionConfig( + qResolution=0.01, + thetaRange=(-12., 12.), + thetaRangeR=(-12., 12.), + fileIdentifier=["610"], + scale=[1], + normalisationFileIdentifier=[], + timeSlize=[300.0] + ) + output_config = options.OutputConfig( + outputFormats=["Rqz.ort"], + outputName='test' + ) + config=options.EOSConfig(self.reader_config, experiment_config, reduction_config, output_config) + reducer = reduction.AmorReduction(config) + reducer.reduce() + + def test_noslicing(self): + experiment_config = options.ExperimentConfig( + chopperPhase=-13.5, + chopperPhaseOffset=-5, + yRange=(11., 41.), + lambdaRange=(2., 15.), + qzRange=(0.005, 0.30), + offSpecular=False, + mu=0, + nu=0, + muOffset=0.0 + ) + reduction_config = options.ReductionConfig( + qResolution=0.01, + thetaRange=(-12., 12.), + thetaRangeR=(-12., 12.), + fileIdentifier=["610", "611", "608,612-613", "609"], + scale=[1], + normalisationFileIdentifier=["614"], + autoscale=(0.005, 0.008) + ) + output_config = options.OutputConfig( + outputFormats=["Rqz.ort"], + outputName='test' + ) + config=options.EOSConfig(self.reader_config, experiment_config, reduction_config, output_config) + reducer = reduction.AmorReduction(config) + reducer.reduce() + # run second time to reuse norm file + reducer = reduction.AmorReduction(config) + reducer.reduce() From 1780916a12b5fc28cba81530dc9be191a0b654d5 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 4 Mar 2024 17:40:00 +0100 Subject: [PATCH 11/17] correct pycache folder name in .gitignore --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index cfe693c..5205740 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,7 @@ *~ *.bak *.log -__py_cache__ +__pycache__ raw .idea test_data From 7274e1bc853772e082d42c4e7d06ea07174b8fef Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 5 Mar 2024 10:19:21 +0100 Subject: [PATCH 12/17] Add profiling to tests and separate read_individual_data into sub-methods --- libeos/file_reader.py | 298 +++++++++++++++++++----------------- tests/test_full_analysis.py | 16 ++ 2 files changed, 173 insertions(+), 141 deletions(-) diff --git a/libeos/file_reader.py b/libeos/file_reader.py index e9f49ca..7adb8e7 100644 --- a/libeos/file_reader.py +++ b/libeos/file_reader.py @@ -123,62 +123,127 @@ class AmorData: #return matr #------------------------------------------------------------------------------------------------- def read_individual_data(self, fileName, norm=False): - pixelLookUp = self.resolve_pixels() - self.hdf = h5py.File(fileName, 'r', swmr=True) if self.readHeaderInfo: - # read general information and first data set - logging.info(f'# meta data from: {self.file_list[0]}') - self.hdf = h5py.File(self.file_list[0], 'r', swmr=True) - - title = self.hdf['entry1/title'][0].decode('utf-8') - proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8') - user_name = self.hdf['entry1/user/name'][0].decode('utf-8') - user_affiliation = 'unknown' - 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 - - # assembling orso header information - self.header.owner = fileio.Person( - name = user_name, - affiliation = user_affiliation, - contact = user_email, - ) - if user_orcid: - self.header.owner.orcid = user_orcid - self.header.experiment = fileio.Experiment( - title = title, - instrument = instrumentName, - start_date = self.start_date, - probe = sourceProbe, - facility = source, - proposalID = proposal_id - ) - self.header.sample = fileio.Sample( - name = sampleName, - model = model, - sample_parameters = None, - ) - self.header.measurement_scheme = 'angle- and energy-dispersive' - - 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 + self.read_header_info() logging.info(f'# data from file: {fileName}') + self.read_individual_header() + + # add header content + 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, + ) + self.header.measurement_instrument_settings.mu = fileio.Value(round(self.mu, 3), 'deg', comment='sample angle to horizon') + self.header.measurement_instrument_settings.nu = fileio.Value(round(self.nu, 3), 'deg', comment='detector angle to horizon') + if norm: + self.header.measurement_additional_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=self.fileDate)) + else: + self.header.measurement_data_files.append(fileio.File(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.kap:6.3f}') + + # TODO: should extract monitor from counts or beam current times time + self.monitor1 = self.ctime + self.monitor2 = self.monitor1 + + self.read_event_stream() + totalNumber = np.shape(self.tof_e)[0] + + self.extract_walltime(norm, totalNumber) + + if True: + self.filter_strange_times() + self.marge_frames() + + self.filter_project_x() + + # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau + # TODO: check for correctness + if not self.config.offSpecular: + self.tof_e -= ( self.delta_e / 180. ) * self.tau + + self.calculate_derived_properties() + + self.filter_qz_range(norm) + + logging.info(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') + + def filter_qz_range(self, norm): + if self.config.qzRange[1]<0.3 and not norm: + self.mask_e = np.logical_and(self.mask_e, + (self.config.qzRange[0]<=self.qz_e) & (self.qz_e<=self.config.qzRange[1])) + self.detZ_e = self.detZ_e[self.mask_e] + self.lamda_e = self.lamda_e[self.mask_e] + self.wallTime_e = self.wallTime_e[self.mask_e] + + def calculate_derived_properties(self): + # lambda + self.lamda_e = 1.e13*self.tof_e*const.hdm/(self.chopperDetectorDistance+self.detXdist_e) + self.lamdaMax = const.lamdaCut+1.e13*self.tau*const.hdm/(self.chopperDetectorDistance+124.) + self.mask_e = np.logical_and(self.mask_e, (self.config.lambdaRange[0]<=self.lamda_e) & ( + self.lamda_e<=self.config.lambdaRange[1])) + # alpha_f + alphaF_e = self.nu-self.mu+self.delta_e + # q_z + if self.config.offSpecular: + 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', + else: + 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' + + def filter_project_x(self): + pixelLookUp = self.resolve_pixels() + # 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 marge_frames(self): + self.tof_e = np.remainder(self.tof_e-self.tofCut+self.tau, self.tau)+self.tofCut # tof shifted to 1 frame + self.tof_e = self.tof_e+self.tau*self.config.chopperPhaseOffset/180. # correction for time offset between chopper pulse and tof zero + + def filter_strange_times(self): + # filter 'strange' tof times > 2 tau + filter_e = (self.tof_e<=2*self.tau) + self.tof_e = self.tof_e[filter_e] + self.pixelID_e = self.pixelID_e[filter_e] + self.wallTime_e = self.wallTime_e[filter_e] + if np.shape(filter_e)[0]-np.shape(self.tof_e)[0]>0.5: + logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]}') + + def extract_walltime(self, norm, totalNumber): + self.wallTime_e = np.empty(totalNumber) + 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] + if not self.startTime and not norm: + self.startTime = self.wallTime_e[0] + self.wallTime_e -= self.startTime + logging.debug(f'wall time from {self.wallTime_e[0]} to {self.wallTime_e[-1]}') + + def read_event_stream(self): + self.tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 + self.pixelID_e = np.array(self.hdf['/entry1/Amor/detector/data/event_id'][:], dtype=int) + self.dataPacket_p = np.array(self.hdf['/entry1/Amor/detector/data/event_index'][:], dtype=np.uint64) + self.dataPacketTime_p = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.uint64)/1e9 + + 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 + try: self.mu = float(np.take(self.hdf['/entry1/Amor/master_parameters/mu/value'], 0)) self.nu = float(np.take(self.hdf['/entry1/Amor/master_parameters/nu/value'], 0)) @@ -216,96 +281,47 @@ class AmorData: # TODO: figure out real stop time.... self.ctime=(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-1] - self.hdf['/entry1/Amor/detector/data/event_time_zero'][0]) / 1.e9 - fileDate = datetime.fromisoformat( self.hdf['/entry1/start_time'][0].decode('utf-8') ) + self.fileDate = datetime.fromisoformat( self.hdf['/entry1/start_time'][0].decode('utf-8') ) - # add header content - 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, + def read_header_info(self): + # read general information and first data set + logging.info(f'# meta data from: {self.file_list[0]}') + self.hdf = h5py.File(self.file_list[0], 'r', swmr=True) + title = self.hdf['entry1/title'][0].decode('utf-8') + proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8') + user_name = self.hdf['entry1/user/name'][0].decode('utf-8') + user_affiliation = 'unknown' + 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 + # assembling orso header information + self.header.owner = fileio.Person( + name=user_name, + affiliation=user_affiliation, + contact=user_email, ) - self.header.measurement_instrument_settings.mu = fileio.Value(round(self.mu, 3), 'deg', comment='sample angle to horizon') - self.header.measurement_instrument_settings.nu = fileio.Value(round(self.nu, 3), 'deg', comment='detector angle to horizon') - if norm: - self.header.measurement_additional_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) - else: - self.header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=fileDate)) - logging.info(f'# mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kap:6.3f}') - - # TODO: should extract monitor from counts or beam current times time - self.monitor1 = self.ctime - self.monitor2 = self.monitor1 - - # read data event streams - tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 - pixelID_e = np.array(self.hdf['/entry1/Amor/detector/data/event_id'][:], dtype=int) - totalNumber = np.shape(tof_e)[0] - - wallTime_e = np.empty(totalNumber) - dataPacket_p = np.array(self.hdf['/entry1/Amor/detector/data/event_index'][:], dtype=np.uint64) - dataPacketTime_p = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.uint64)/1e9 - #for i, index in enumerate(dataPacket_p): - # wallTime_e[index:] = dataPacketTime_p[i] - for i in range(len(dataPacket_p)-1): - wallTime_e[dataPacket_p[i]:dataPacket_p[i+1]] = dataPacketTime_p[i] - wallTime_e[dataPacket_p[-1]:] = dataPacketTime_p[-1] - if not self.startTime and not norm: - self.startTime = wallTime_e[0] - wallTime_e -= self.startTime - logging.debug(f'wall time from {wallTime_e[0]} to {wallTime_e[-1]}') - - # filter 'strange' tof times > 2 tau - if True: - filter_e = (tof_e <= 2*self.tau) - tof_e = tof_e[filter_e] - pixelID_e = pixelID_e[filter_e] - wallTime_e = wallTime_e[filter_e] - if np.shape(filter_e)[0]-np.shape(tof_e)[0] > 0.5 : - logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(tof_e)[0]}') - tof_e = np.remainder( tof_e - self.tofCut + self.tau, self.tau) + self.tofCut # tof shifted to 1 frame - tof_e = tof_e + self.tau * self.config.chopperPhaseOffset / 180. # correction for time offset between chopper pulse and tof zero - - # resolve pixel ID into y and z indicees, x position and angle - (detY_e, detZ_e, detXdist_e, delta_e) = pixelLookUp[np.int_(pixelID_e)-1,:].T - - # define mask and filter y range - mask_e = (self.config.yRange[0] <= detY_e) & (detY_e <= self.config.yRange[1]) - - # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau - # TODO: check for correctness - if not self.config.offSpecular: - tof_e -= ( delta_e / 180. ) * self.tau - - # lambda - lamda_e = 1.e13 * tof_e * const.hdm / (self.chopperDetectorDistance + detXdist_e) - self.lamdaMax = const.lamdaCut + 1.e13 * self.tau * const.hdm / (self.chopperDetectorDistance + 124.) - mask_e = np.logical_and(mask_e, (self.config.lambdaRange[0] <= lamda_e) & (lamda_e <= self.config.lambdaRange[1])) - - # alpha_f - alphaF_e = self.nu - self.mu + delta_e - - # q_z - if self.config.offSpecular: - alphaI = self.kap + self.kad + self.mu - qz_e = 2 * np.pi * ( np.sin( np.deg2rad(alphaF_e) ) + np.sin( np.deg2rad( alphaI ) ) ) / lamda_e - qx_e = 2 * np.pi * ( np.cos( np.deg2rad(alphaF_e) ) - np.cos( np.deg2rad( alphaI ) ) ) / lamda_e - self.header.measurement_scheme = 'energy-dispersive', - else: - qz_e = 4 * np.pi * np.sin( np.deg2rad(alphaF_e) ) / lamda_e - # qx_e = 0. - self.header.measurement_scheme = 'angle- and energy-dispersive' - - # filter q_z range - if self.config.qzRange[1] < 0.3 and not norm: - mask_e = np.logical_and(mask_e, (self.config.qzRange[0] <= qz_e) & (qz_e <= self.config.qzRange[1])) - - self.detZ_e = detZ_e[mask_e] - self.lamda_e = lamda_e[mask_e] - self.wallTime_e = wallTime_e[mask_e] - - logging.info(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') + if user_orcid: + self.header.owner.orcid = user_orcid + self.header.experiment = fileio.Experiment( + title=title, + instrument=instrumentName, + start_date=self.start_date, + probe=sourceProbe, + facility=source, + proposalID=proposal_id + ) + self.header.sample = fileio.Sample( + name=sampleName, + model=model, + sample_parameters=None, + ) + self.header.measurement_scheme = 'angle- and energy-dispersive' diff --git a/tests/test_full_analysis.py b/tests/test_full_analysis.py index 29606c9..44a6878 100644 --- a/tests/test_full_analysis.py +++ b/tests/test_full_analysis.py @@ -1,4 +1,5 @@ import os +import cProfile from unittest import TestCase from libeos import options, reduction, logconfig @@ -6,12 +7,22 @@ logconfig.setup_logging() logconfig.update_loglevel(True, False) class FullAmorTest(TestCase): + @classmethod + def setUpClass(cls): + cls.pr = cProfile.Profile() + + @classmethod + def tearDownClass(cls): + cls.pr.dump_stats("profile_test.prof") + def setUp(self): + self.pr.enable() self.reader_config = options.ReaderConfig( year=2023, dataPath=os.path.join('..', "test_data")) def tearDown(self): + self.pr.disable() for fi in ['test.Rqz.ort', '614.norm']: try: os.unlink(os.path.join(self.reader_config.dataPath, fi)) @@ -46,6 +57,11 @@ class FullAmorTest(TestCase): outputName='test' ) config=options.EOSConfig(self.reader_config, experiment_config, reduction_config, output_config) + # run three times to get similar timing to noslicing runs + reducer = reduction.AmorReduction(config) + reducer.reduce() + reducer = reduction.AmorReduction(config) + reducer.reduce() reducer = reduction.AmorReduction(config) reducer.reduce() From 736ce0b2e92694ba6ec12ade7116e1f2a81af8a6 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 5 Mar 2024 14:09:09 +0100 Subject: [PATCH 13/17] Add alternative histogram function, but did not improve performance for the size of lambda binning used --- libeos/reduction.py | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/libeos/reduction.py b/libeos/reduction.py index 44c6073..ad2f34c 100644 --- a/libeos/reduction.py +++ b/libeos/reduction.py @@ -385,3 +385,46 @@ class AmorReduction: return qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz + + @staticmethod + def histogram2d_lz(lamda_e, detZ_e, bins): + """ + Perform binning operation equivalent to numpy bin2d for the sepcial case + of the second dimension using integer positions (pre-defined pixels). + Based on the devide_bin algorithm below. + """ + dimension = bins[1].shape[0]-1 + if not (np.array(bins[1])==np.arange(dimension+1)).all(): + raise ValueError("histogram2d_lz requires second bin dimension to be contigous integer range") + binning = AmorReduction.devide_bin(lamda_e, detZ_e.astype(np.int64), bins[0], dimension) + return np.array(binning), bins[0], bins[1] + + @staticmethod + def devide_bin(lambda_e, position_e, lamda_edges, dimension): + ''' + Use a divide and conquer strategy to bin the data. For the actual binning the + numpy bincount function is used, as it is much faster than histogram for + counting of integer values. + + :param lambda_e: Array of wavelength for each event + :param position_e: Array of positional indices for each event + :param lamda_edges: The edges of bins to be used for the histogram + :param dimension: position number of buckets in output arrray + + :return: 2D list of dimensions (lambda, x) of counts + ''' + if len(lambda_e)==0: + # no more events in range, return empty bins + return [np.zeros(dimension, dtype=int).tolist()]*(len(lamda_edges)-1) + if len(lamda_edges)==2: + # deepest recursion reached, all items should be within the two ToF edges + return [np.bincount(position_e, minlength=dimension).tolist()] + # split all events into two time of flight regions + split_idx = len(lamda_edges)//2 + left_region = lambda_e Date: Tue, 5 Mar 2024 14:14:51 +0100 Subject: [PATCH 14/17] Speed up merge_frames by 45% using combined scalar offsets instead of full array operations --- libeos/file_reader.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/libeos/file_reader.py b/libeos/file_reader.py index 7adb8e7..64de609 100644 --- a/libeos/file_reader.py +++ b/libeos/file_reader.py @@ -160,7 +160,7 @@ class AmorData: if True: self.filter_strange_times() - self.marge_frames() + self.merge_frames() self.filter_project_x() @@ -209,9 +209,9 @@ class AmorData: # define mask and filter y range self.mask_e = (self.config.yRange[0]<=detY_e) & (detY_e<=self.config.yRange[1]) - def marge_frames(self): - self.tof_e = np.remainder(self.tof_e-self.tofCut+self.tau, self.tau)+self.tofCut # tof shifted to 1 frame - self.tof_e = self.tof_e+self.tau*self.config.chopperPhaseOffset/180. # correction for time offset between chopper pulse and tof zero + def merge_frames(self): + total_offset = self.tofCut+self.tau*self.config.chopperPhaseOffset/180. + self.tof_e = np.remainder(self.tof_e-(self.tofCut-self.tau), self.tau)+total_offset # tof shifted to 1 frame def filter_strange_times(self): # filter 'strange' tof times > 2 tau From 43313841d6723a2444319b17ed91cd89f38e8d5a Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 5 Mar 2024 14:33:41 +0100 Subject: [PATCH 15/17] minor --- libeos/file_reader.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/libeos/file_reader.py b/libeos/file_reader.py index 64de609..87fc1e2 100644 --- a/libeos/file_reader.py +++ b/libeos/file_reader.py @@ -156,7 +156,7 @@ class AmorData: self.read_event_stream() totalNumber = np.shape(self.tof_e)[0] - self.extract_walltime(norm, totalNumber) + self.extract_walltime(norm) if True: self.filter_strange_times() @@ -185,7 +185,7 @@ class AmorData: def calculate_derived_properties(self): # lambda - self.lamda_e = 1.e13*self.tof_e*const.hdm/(self.chopperDetectorDistance+self.detXdist_e) + self.lamda_e = (1.e13*self.tof_e*const.hdm)/(self.chopperDetectorDistance+self.detXdist_e) self.lamdaMax = const.lamdaCut+1.e13*self.tau*const.hdm/(self.chopperDetectorDistance+124.) self.mask_e = np.logical_and(self.mask_e, (self.config.lambdaRange[0]<=self.lamda_e) & ( self.lamda_e<=self.config.lambdaRange[1])) @@ -222,7 +222,8 @@ class AmorData: if np.shape(filter_e)[0]-np.shape(self.tof_e)[0]>0.5: logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]}') - def extract_walltime(self, norm, totalNumber): + def extract_walltime(self, norm): + totalNumber = np.shape(self.tof_e)[0] self.wallTime_e = np.empty(totalNumber) 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] From 146743673ca60df170440ef38dc95dd2b113e421 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 5 Mar 2024 15:49:15 +0100 Subject: [PATCH 16/17] Add numba implementation of some event functions, yielding x2 speed improvement, also some small general speedups --- libeos/file_reader.py | 53 +++++++++++++++++++++++++---------- libeos/nb_helpers.py | 65 +++++++++++++++++++++++++++++++++++++++++++ requirements.txt | 1 + 3 files changed, 104 insertions(+), 15 deletions(-) create mode 100644 libeos/nb_helpers.py diff --git a/libeos/file_reader.py b/libeos/file_reader.py index 87fc1e2..9824c9d 100644 --- a/libeos/file_reader.py +++ b/libeos/file_reader.py @@ -14,6 +14,11 @@ from .header import Header from .instrument import Detector from .options import ExperimentConfig, ReaderConfig +try: + from . import nb_helpers +except Exception: + nb_helpers = None + class AmorData: """read meta-data and event streams from .hdf file(s), apply filters and conversions""" @@ -184,9 +189,16 @@ class AmorData: self.wallTime_e = self.wallTime_e[self.mask_e] def calculate_derived_properties(self): - # lambda - self.lamda_e = (1.e13*self.tof_e*const.hdm)/(self.chopperDetectorDistance+self.detXdist_e) self.lamdaMax = const.lamdaCut+1.e13*self.tau*const.hdm/(self.chopperDetectorDistance+124.) + if nb_helpers and not self.config.offSpecular: + 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 + ) + return + # 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])) # alpha_f @@ -194,24 +206,32 @@ class AmorData: # q_z if self.config.offSpecular: 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.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', else: - self.qz_e = 4*np.pi*np.sin(np.deg2rad(alphaF_e))/self.lamda_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' def filter_project_x(self): pixelLookUp = self.resolve_pixels() - # 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]) + 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 merge_frames(self): total_offset = self.tofCut+self.tau*self.config.chopperPhaseOffset/180. - self.tof_e = np.remainder(self.tof_e-(self.tofCut-self.tau), self.tau)+total_offset # tof shifted to 1 frame + 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_strange_times(self): # filter 'strange' tof times > 2 tau @@ -223,11 +243,14 @@ class AmorData: logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]}') def extract_walltime(self, norm): - totalNumber = np.shape(self.tof_e)[0] - self.wallTime_e = np.empty(totalNumber) - 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] + if nb_helpers: + self.wallTime_e = nb_helpers.extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p) + else: + totalNumber = np.shape(self.tof_e)[0] + self.wallTime_e = np.empty(totalNumber) + 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] if not self.startTime and not norm: self.startTime = self.wallTime_e[0] self.wallTime_e -= self.startTime diff --git a/libeos/nb_helpers.py b/libeos/nb_helpers.py new file mode 100644 index 0000000..8f69824 --- /dev/null +++ b/libeos/nb_helpers.py @@ -0,0 +1,65 @@ +import numba as nb +import numpy as np + +@nb.jit(nb.float64[:](nb.float64[:], nb.float64, nb.float64, nb.float64), + nopython=True, parallel=True, cache=True) +def merge_frames(tof_e, tofCut, tau, total_offset): + # fast implementation of the merging function used in file_reader.AmorData.merge_frames + tof_e_out = np.empty(tof_e.shape, dtype=np.float64) + dt = (tofCut-tau) + for ti in nb.prange(tof_e.shape[0]): + tof_e_out[ti] = ((tof_e[ti]-dt)%tau)+total_offset # tof shifted to 1 frame + return tof_e + +@nb.jit(nb.float64[:](nb.float64[:], nb.uint64[:], nb.float64[:]), + nopython=True, parallel=True, cache=True) +def extract_walltime(tof_e, dataPacket_p, dataPacketTime_p): + # assigning every event the wall time of the event packet (absolute time of pulse ?start?) + totalNumber = np.shape(tof_e)[0] + wallTime_e = np.empty(totalNumber, dtype=np.float64) + for i in nb.prange(len(dataPacket_p)-1): + for j in range(dataPacket_p[i], dataPacket_p[i+1]): + wallTime_e[j] = dataPacketTime_p[i] + for j in range(dataPacket_p[-1], totalNumber): + wallTime_e[j] = dataPacketTime_p[-1] + return wallTime_e + +@nb.jit(nb.types.Tuple((nb.int64[:], nb.float64[:], nb.float64[:], nb.boolean[:])) + (nb.float64[:, :], nb.int64[:], nb.int64, nb.int64), + nopython=True, parallel=True, cache=True) +def filter_project_x(pixelLookUp, pixelID_e, ymin, ymax): + # project events on z-axis and create filter for events outside of y-range + events = pixelID_e.shape[0] + detY_e = np.empty(events, dtype=np.int64) + detZ_e = np.empty(events, dtype=np.int64) + detXdist_e = np.empty(events, dtype=np.float64) + delta_e = np.empty(events, dtype=np.float64) + mask_e = np.empty(events, dtype=nb.boolean) + for i in nb.prange(events): + # resolve pixel ID into y and z indicees, x position and angle + detY_e[i] = pixelLookUp[pixelID_e[i]-1, 0] + detZ_e[i] = pixelLookUp[pixelID_e[i]-1, 1] + detXdist_e[i] = pixelLookUp[pixelID_e[i]-1, 2] + delta_e[i] = pixelLookUp[pixelID_e[i]-1, 3] + # define mask and filter y range + mask_e[i] = (ymin<=detY_e[i]) & (detY_e[i]<=ymax) + return (detZ_e, detXdist_e, delta_e, mask_e) + +@nb.jit(nb.types.Tuple((nb.float64[:], nb.float64[:], nb.boolean[:])) + (nb.float64[:], nb.float64[:], nb.float64[:], nb.boolean[:], + nb.float64, nb.float64, nb.float64, nb.float64, nb.float64, nb.float64), + nopython=True, parallel=True, cache=True) +def calculate_derived_properties_focussing(tof_e, detXdist_e, delta_e, mask_e, + lmin, lmax, nu, mu, chopperDetectorDistance, hdm): + events = tof_e.shape[0] + alphaF_e = np.empty(events, dtype=np.float64) + lamda_e = np.empty(events, dtype=np.float64) + qz_e = np.empty(events, dtype=np.float64) + mask_e_out = np.empty(events, dtype=nb.boolean) + denom_f1 = 1.e13*hdm + for i in nb.prange(events): + lamda_e[i] = denom_f1*tof_e[i]/(chopperDetectorDistance+detXdist_e[i]) + mask_e_out[i] = mask_e[i] & ((lmin<=lamda_e[i]) & (lamda_e[i]<=lmax)) + alphaF_e[i] = nu-mu+delta_e[i] + qz_e[i] = 4*np.pi*np.sin(np.deg2rad(alphaF_e[i]))/lamda_e[i] + return lamda_e, qz_e, mask_e \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index c92ee06..b86897a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ numpy h5py orsopy +numba From 11c24a44e12ebdf05d21194c99185770fd42f5f8 Mon Sep 17 00:00:00 2001 From: Jochen Stahn <57442805+jochenstahn@users.noreply.github.com> Date: Wed, 6 Mar 2024 07:51:44 +0100 Subject: [PATCH 17/17] added Artur as author --- neos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/neos.py b/neos.py index 8264173..61a2232 100644 --- a/neos.py +++ b/neos.py @@ -2,7 +2,8 @@ """ eos reduces measurements performed on Amor@SINQ, PSI -Author: Jochen Stahn +Author: Jochen Stahn (algorithms, python draft), + Artur Glavic (structuring and optimisation of code) conventions (not strictly followed, yet): - array names end with the suffix '_x[y]' with the meaning @@ -13,15 +14,14 @@ conventions (not strictly followed, yet): _z = detector z _lz = (lambda, detector z) _q = q_z -- to come """ + import logging from libeos.command_line import command_line_options from libeos.logconfig import setup_logging from libeos.reduction import AmorReduction - #===================================================================================================== # TODO: # - calculate resolution using the chopperPhase