Merge pull request #2 from jochenstahn/input_artur
Refactor, first tests and speedup
This commit is contained in:
Jochen Stahn
+4
-2
@@ -1,5 +1,7 @@
|
||||
*~
|
||||
*.bak
|
||||
__py_cache__
|
||||
*.log
|
||||
__pycache__
|
||||
raw
|
||||
|
||||
.idea
|
||||
test_data
|
||||
|
||||
@@ -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'
|
||||
@@ -0,0 +1,193 @@
|
||||
import argparse
|
||||
from datetime import datetime
|
||||
|
||||
from .logconfig import update_loglevel
|
||||
from .options import ReaderConfig, EOSConfig, ExperimentConfig, OutputConfig, ReductionConfig
|
||||
|
||||
|
||||
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 <interval> ,[<start> [,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)
|
||||
|
||||
def command_line_options():
|
||||
clas = commandLineArgs()
|
||||
update_loglevel(clas.verbose, clas.debug)
|
||||
|
||||
reader_config = ReaderConfig(
|
||||
year=clas.year,
|
||||
dataPath=clas.dataPath)
|
||||
experiment_config = ExperimentConfig(
|
||||
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,
|
||||
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 EOSConfig(reader_config, experiment_config, reduction_config, output_config)
|
||||
@@ -0,0 +1,6 @@
|
||||
"""
|
||||
Constants used in data reduction.
|
||||
"""
|
||||
|
||||
hdm = 6.626176e-34/1.674928e-27 # h / m
|
||||
lamdaCut = 2.5 # Aa
|
||||
@@ -0,0 +1,351 @@
|
||||
import logging
|
||||
import os
|
||||
import subprocess
|
||||
import sys
|
||||
from datetime import datetime
|
||||
from typing import List
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
from orsopy import fileio
|
||||
|
||||
from . import const
|
||||
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"""
|
||||
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:str, norm=False):
|
||||
self.startTime = reader_config.startTime
|
||||
self.header = header
|
||||
self.config = config
|
||||
self.reader_config = reader_config
|
||||
self.expand_file_list(short_notation)
|
||||
self.read_data(norm=norm)
|
||||
|
||||
#-------------------------------------------------------------------------------------------------
|
||||
def read_data(self, norm=False):
|
||||
self.file_list = []
|
||||
for number in self.data_file_numbers:
|
||||
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.reader_config.year}n{number:06d}.hdf'
|
||||
if os.path.exists(os.path.join(self.reader_config.dataPath,fileName)):
|
||||
path = self.reader_config.dataPath
|
||||
elif os.path.exists(fileName):
|
||||
path = '.'
|
||||
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 = 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 os.path.join(path, fileName)
|
||||
#-------------------------------------------------------------------------------------------------
|
||||
def expand_file_list(self, short_notation):
|
||||
"""Evaluate string entry for file number lists"""
|
||||
#log().debug('Executing get_flist')
|
||||
file_list=[]
|
||||
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)]
|
||||
self.data_file_numbers=sorted(file_list)
|
||||
#-------------------------------------------------------------------------------------------------
|
||||
def resolve_pixels(self):
|
||||
"""determine spatial coordinats and angles from pixel number"""
|
||||
nPixel = Detector.nWires * Detector.nStripes * Detector.nBlades
|
||||
pixelID = np.arange(nPixel)
|
||||
(bladeNr, bPixel) = np.divmod(pixelID, Detector.nWires * Detector.nStripes)
|
||||
(bZi, detYi) = np.divmod(bPixel, Detector.nStripes) # z index on blade, y index on detector
|
||||
detZi = bladeNr * Detector.nWires + bZi # z index on detector
|
||||
detX = bZi * Detector.dX # x position in detector
|
||||
# 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
|
||||
#-------------------------------------------------------------------------------------------------
|
||||
def read_individual_data(self, fileName, norm=False):
|
||||
self.hdf = h5py.File(fileName, 'r', swmr=True)
|
||||
|
||||
if self.readHeaderInfo:
|
||||
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)
|
||||
|
||||
if True:
|
||||
self.filter_strange_times()
|
||||
self.merge_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):
|
||||
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
|
||||
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()
|
||||
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.
|
||||
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
|
||||
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):
|
||||
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
|
||||
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))
|
||||
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
|
||||
self.fileDate = datetime.fromisoformat( self.hdf['/entry1/start_time'][0].decode('utf-8') )
|
||||
|
||||
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,
|
||||
)
|
||||
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'
|
||||
|
||||
@@ -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 . 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)
|
||||
@@ -0,0 +1,68 @@
|
||||
"""
|
||||
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.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(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*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(Detector.nBlades-1):
|
||||
delta_grid = np.concatenate((delta_grid, blade_grid), axis=None)
|
||||
blade_grid = blade_grid + bladeAngle
|
||||
delta_grid = delta_grid[delta_grid<blade_grid[0]-0.5*stepWidth]
|
||||
delta_grid = np.concatenate((delta_grid, blade_grid), axis=None)
|
||||
|
||||
return -np.flip(delta_grid) + 0.5*Detector.nBlades * bladeAngle
|
||||
@@ -0,0 +1,43 @@
|
||||
"""
|
||||
Setup for the logging of eos.
|
||||
"""
|
||||
import sys
|
||||
import logging
|
||||
import logging.handlers
|
||||
|
||||
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 update_loglevel(verbose=False, debug=False):
|
||||
if verbose:
|
||||
logging.getLogger().handlers[0].setLevel(logging.INFO)
|
||||
if debug:
|
||||
console = logging.getLogger().handlers[0]
|
||||
console.setLevel(logging.DEBUG)
|
||||
formatter = logging.Formatter('%(levelname).1s %(message)s')
|
||||
console.setFormatter(formatter)
|
||||
@@ -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
|
||||
@@ -0,0 +1,51 @@
|
||||
"""
|
||||
Classes for stroing various configurations needed for reduction.
|
||||
"""
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Optional, Tuple
|
||||
|
||||
@dataclass
|
||||
class ReaderConfig:
|
||||
year: int
|
||||
dataPath: str
|
||||
startTime: Optional[float] = 0
|
||||
|
||||
@dataclass
|
||||
class ExperimentConfig:
|
||||
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
|
||||
muOffset: Optional[float] = None
|
||||
offSpecular: bool = False
|
||||
|
||||
@dataclass
|
||||
class ReductionConfig:
|
||||
qResolution: float
|
||||
thetaRange: Tuple[float, float]
|
||||
thetaRangeR: 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
|
||||
|
||||
@dataclass
|
||||
class OutputConfig:
|
||||
outputFormats: list
|
||||
outputName: str
|
||||
|
||||
@dataclass
|
||||
class EOSConfig:
|
||||
reader: ReaderConfig
|
||||
experiment: ExperimentConfig
|
||||
reductoin: ReductionConfig
|
||||
output: OutputConfig
|
||||
@@ -0,0 +1,430 @@
|
||||
import logging
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
from orsopy import fileio
|
||||
|
||||
from .command_line import expand_file_list
|
||||
from .file_reader import AmorData
|
||||
from .header import Header
|
||||
from .options import EOSConfig
|
||||
from .instrument import Grid
|
||||
|
||||
class AmorReduction:
|
||||
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()
|
||||
|
||||
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}')
|
||||
|
||||
# load or create normalisation matrix
|
||||
if self.reduction_config.normalisationFileIdentifier:
|
||||
self.create_normalisation_map(self.reduction_config.normalisationFileIdentifier[0])
|
||||
else:
|
||||
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:
|
||||
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:
|
||||
self.subtract = False
|
||||
|
||||
# load measurement data and do the reduction
|
||||
self.datasetsRqz = []
|
||||
self.datasetsRlt = []
|
||||
for i, short_notation in enumerate(self.reduction_config.fileIdentifier):
|
||||
self.read_file_block(i, short_notation)
|
||||
|
||||
# output
|
||||
logging.warning('output:')
|
||||
|
||||
if 'Rqz.ort' in self.output_config.outputFormats:
|
||||
self.save_Rqz()
|
||||
|
||||
if 'Rlt.ort' in self.output_config.outputFormats:
|
||||
self.save_Rtl()
|
||||
|
||||
def read_file_block(self, i, short_notation):
|
||||
logging.warning('reading input:')
|
||||
self.header.measurement_data_files = []
|
||||
self.file_reader = AmorData(header=self.header,
|
||||
reader_config=self.reader_config,
|
||||
config=self.experiment_config,
|
||||
short_notation=short_notation)
|
||||
if self.reduction_config.timeSlize:
|
||||
self.read_timeslices(i)
|
||||
else:
|
||||
self.read_unsliced(i)
|
||||
|
||||
def read_unsliced(self, i):
|
||||
lamda_e = self.file_reader.lamda_e
|
||||
detZ_e = self.file_reader.detZ_e
|
||||
qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, self.mask_lz = self.project_on_lz(
|
||||
self.file_reader, self.norm_lz, self.normAngle, lamda_e, detZ_e)
|
||||
try:
|
||||
ref_lz *= self.reduction_config.scale[i]
|
||||
err_lz *= self.reduction_config.scale[i]
|
||||
except IndexError:
|
||||
ref_lz *= self.reduction_config.scale[-1]
|
||||
err_lz *= self.reduction_config.scale[-1]
|
||||
if 'Rqz.ort' in self.output_config.outputFormats:
|
||||
headerRqz = self.header.orso_header()
|
||||
headerRqz.data_set = f'Nr {i} : mu = {self.file_reader.mu:6.3f} deg'
|
||||
|
||||
# projection on q-grid
|
||||
q_q, R_q, dR_q, dq_q = self.project_on_qz(qz_lz, ref_lz, err_lz, res_lz, self.norm_lz, self.mask_lz)
|
||||
|
||||
filter_q = np.where((self.experiment_config.qzRange[0]<q_q) & (q_q<self.experiment_config.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 self.reduction_config.autoscale:
|
||||
if i==0:
|
||||
R_q, dR_q = self.autoscale(q_q, R_q, dR_q)
|
||||
else:
|
||||
pRq_z = self.datasetsRqz[i-1].data[:, 1]
|
||||
pdRq_z = self.datasetsRqz[i-1].data[:, 2]
|
||||
R_q, dR_q = self.autoscale(q_q, R_q, dR_q, pRq_z, pdRq_z)
|
||||
|
||||
if self.subtract:
|
||||
if len(q_q)==len(self.sq_q):
|
||||
R_q -= self.sR_q
|
||||
dR_q = np.sqrt(dR_q**2+self.sdR_q**2)
|
||||
else:
|
||||
logging.warning(
|
||||
f'backgroung file {self.sFileName} not compatible with q_z scale ({len(self.sq_q)} vs. {len(q_q)})')
|
||||
|
||||
data = np.array([q_q, R_q, dR_q, dq_q]).T
|
||||
orso_data = fileio.OrsoDataset(headerRqz, data)
|
||||
self.datasetsRqz.append(orso_data)
|
||||
if 'Rlt.ort' in self.output_config.outputFormats:
|
||||
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 = file_reader.data_source
|
||||
|
||||
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,
|
||||
self.norm_lz.T,
|
||||
np.where(self.mask_lz, 1, 0).T
|
||||
):
|
||||
data = np.array(list(item)).T
|
||||
headerRlt = self.header.orso_header(columns=columns)
|
||||
headerRlt.data_set = f'dataset_{i}_{j+1} : alpha_f = {theta_lz[0, j]:6.3f} deg'
|
||||
orso_data = fileio.OrsoDataset(headerRlt, data)
|
||||
self.datasetsRlt.append(orso_data)
|
||||
j += 1
|
||||
|
||||
def read_timeslices(self, i):
|
||||
wallTime_e = self.file_reader.wallTime_e
|
||||
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 ti, time in enumerate(np.arange(start, stop, interval)):
|
||||
logging.warning(f' time slize {ti:4d}')
|
||||
|
||||
filter_e = np.where((time<wallTime_e) & (wallTime_e<time+interval), True, False)
|
||||
lamda_e = self.file_reader.lamda_e[filter_e]
|
||||
detZ_e = self.file_reader.detZ_e[filter_e]
|
||||
|
||||
qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz = self.project_on_lz(
|
||||
self.file_reader, self.norm_lz, self.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, self.norm_lz, mask_lz)
|
||||
|
||||
filter_q = np.where((self.experiment_config.qzRange[0]<q_q) & (q_q<self.experiment_config.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 self.reduction_config.autoscale:
|
||||
R_q, dR_q = self.autoscale(q_q, R_q, dR_q)
|
||||
|
||||
if self.subtract:
|
||||
if len(q_q)==len(self.sq_q):
|
||||
R_q -= self.sR_q
|
||||
dR_q = np.sqrt(dR_q**2+self.sdR_q**2)
|
||||
else:
|
||||
self.subtract = False
|
||||
logging.warning(
|
||||
f'background file {self.sFileName} not compatible with q_z scale ({len(self.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 = self.header.orso_header(
|
||||
extra_columns=[fileio.Column('time', 's', 'time relative to start of measurement series')])
|
||||
headerRqz.data_set = f'{i}_{ti}: time = {time:8.1f} s to {time+interval:8.1f} s'
|
||||
orso_data = fileio.OrsoDataset(headerRqz, data)
|
||||
self.datasetsRqz.append(orso_data)
|
||||
# reset normal logging behavior
|
||||
logging.StreamHandler.terminator = "\n"
|
||||
logging.warning(f' time slize {ti:4d}, done')
|
||||
|
||||
def save_Rqz(self):
|
||||
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, fname, data_separator='\n', comment=theSecondLine)
|
||||
|
||||
def save_Rtl(self):
|
||||
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, fname, data_separator='\n', comment=theSecondLine)
|
||||
|
||||
def autoscale(self, q_q, R_q, dR_q, pR_q=[], pdR_q=[]):
|
||||
autoscale = self.reduction_config.autoscale
|
||||
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
|
||||
|
||||
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):
|
||||
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 \'{fname}\' 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 create_normalisation_map(self, short_notation):
|
||||
dataPath = self.reader_config.dataPath
|
||||
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]}'
|
||||
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')
|
||||
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
|
||||
else:
|
||||
logging.info(f'# normalisation matrix: using the files {normalisation_list}')
|
||||
fromHDF = AmorData(header=self.header,
|
||||
reader_config=self.reader_config,
|
||||
config=self.experiment_config,
|
||||
short_notation=short_notation, norm=True)
|
||||
self.normAngle = fromHDF.nu - fromHDF.mu
|
||||
lamda_e = fromHDF.lamda_e
|
||||
detZ_e = fromHDF.detZ_e
|
||||
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 = 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)
|
||||
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 = {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', 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
|
||||
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.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 ) )
|
||||
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
|
||||
|
||||
|
||||
@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<lamda_edges[split_idx]
|
||||
left_list = AmorReduction.devide_bin(lambda_e[left_region], position_e[left_region],
|
||||
lamda_edges[:split_idx+1], dimension)
|
||||
right_region = np.logical_not(left_region)
|
||||
right_list = AmorReduction.devide_bin(lambda_e[right_region], position_e[right_region],
|
||||
lamda_edges[split_idx:], dimension)
|
||||
return left_list+right_list
|
||||
@@ -1,9 +1,9 @@
|
||||
#!/usr/bin/env python
|
||||
#-*- coding: utf-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
|
||||
@@ -14,909 +14,33 @@ conventions (not strictly followed, yet):
|
||||
_z = detector z
|
||||
_lz = (lambda, detector z)
|
||||
_q = q_z
|
||||
- to come
|
||||
"""
|
||||
|
||||
__version__ = '2.0'
|
||||
__date__ = '2024-03-01'
|
||||
|
||||
import os
|
||||
import sys
|
||||
import subprocess
|
||||
import logging
|
||||
import argparse
|
||||
import numpy as np
|
||||
from datetime import datetime
|
||||
import h5py
|
||||
from orsopy import fileio
|
||||
import platform
|
||||
|
||||
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
|
||||
# - 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 <interval> ,[<start> [,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")
|
||||
|
||||
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
|
||||
print(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
|
||||
|
||||
print(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))
|
||||
print(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
|
||||
#print(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 :
|
||||
print(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]
|
||||
|
||||
print(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'):
|
||||
print(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:
|
||||
print(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:
|
||||
print(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:
|
||||
print(f'# automatic scaling not possible')
|
||||
scale = 1.
|
||||
R_q /= scale
|
||||
dR_q /= scale
|
||||
#print(f'# scaling factor = {scale}')
|
||||
|
||||
return R_q, dR_q
|
||||
#=====================================================================================================
|
||||
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
|
||||
#=====================================================================================================
|
||||
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:
|
||||
print(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<blade_grid[0]-0.5*stepWidth]
|
||||
delta_grid = np.concatenate((delta_grid, blade_grid), axis=None)
|
||||
|
||||
return -np.flip(delta_grid) + 0.5*self.det.nBlades * bladeAngle
|
||||
|
||||
#=====================================================================================================
|
||||
def main():
|
||||
global startTime, grid, clas, header
|
||||
clas = commandLineArgs()
|
||||
grid = Grid()
|
||||
header = Header()
|
||||
startTime = 0
|
||||
if not os.path.exists(f'{clas.dataPath}'):
|
||||
os.system(f'mkdir {clas.dataPath}')
|
||||
fromHDF = AmorData()
|
||||
|
||||
print('\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.reduction.corrections.append('normalisation with \'additional files\'')
|
||||
else:
|
||||
normalise = False
|
||||
norm_lz = grid.lz()
|
||||
normAngle = 1.
|
||||
setup_logging()
|
||||
logging.warning('######## eos - data reduction for Amor ########')
|
||||
|
||||
print('# normalisation matrix: none requested')
|
||||
# read command line arguments and generate classes holding configuration parameters
|
||||
config = command_line_options()
|
||||
# Create reducer with these arguments
|
||||
reducer = AmorReduction(config)
|
||||
# Perform actual reduction
|
||||
reducer.reduce()
|
||||
|
||||
# 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}')
|
||||
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):
|
||||
print('# reading input:')
|
||||
header.measurement_data_files = []
|
||||
fromHDF.read_data(short_notation)
|
||||
logging.info('######## eos - finished ########')
|
||||
|
||||
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]
|
||||
for i, time in enumerate(np.arange(start, stop, interval)):
|
||||
print(f'# time slize {i:4d}', end='\r')
|
||||
|
||||
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)
|
||||
q_q, R_q, dR_q, dq_q = project_on_qz(qz_lz, ref_lz, err_lz, res_lz, norm_lz, mask_lz)
|
||||
|
||||
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)
|
||||
|
||||
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
|
||||
print(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)
|
||||
print('')
|
||||
|
||||
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)
|
||||
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)
|
||||
|
||||
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:
|
||||
if i == 0:
|
||||
R_q, dR_q = autoscale(q_q, R_q, dR_q)
|
||||
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, 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:
|
||||
print(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
|
||||
print('# output:')
|
||||
|
||||
if 'Rqz.ort' in output_format_list(clas.outputFormat):
|
||||
print(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')
|
||||
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()
|
||||
|
||||
@@ -0,0 +1,4 @@
|
||||
numpy
|
||||
h5py
|
||||
orsopy
|
||||
numba
|
||||
@@ -0,0 +1,98 @@
|
||||
import os
|
||||
import cProfile
|
||||
from unittest import TestCase
|
||||
from libeos import options, reduction, logconfig
|
||||
|
||||
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))
|
||||
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)
|
||||
# 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()
|
||||
|
||||
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()
|
||||
Reference in New Issue
Block a user