Compare commits
18 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
01873d60db | ||
| 654251a194 | |||
| 9e30970d9b | |||
| 7b6f2045cc | |||
| e3d15a1142 | |||
|
|
98244327eb | ||
| 61ba1cf084 | |||
| 54f2169888 | |||
| 24cc0a4287 | |||
| a3d342e5d7 | |||
| c92e4e4d17 | |||
| 4d74237669 | |||
| 6144200551 | |||
| f2c681ffba | |||
| a56f29191d | |||
| 6597ca22d1 | |||
| ce31dcb190 | |||
| 41c5218413 |
@@ -2,5 +2,5 @@
|
||||
Package to handle data redction at AMOR instrument to be used by __main__.py script.
|
||||
"""
|
||||
|
||||
__version__ = '3.0.6'
|
||||
__date__ = '2025-11-05'
|
||||
__version__ = '3.1.0'
|
||||
__date__ = '2026-02-20'
|
||||
|
||||
@@ -29,7 +29,7 @@ except ImportError:
|
||||
AMOR_LOCAL_TIMEZONE = zoneinfo.ZoneInfo(key='Europe/Zurich')
|
||||
|
||||
if platform.node().startswith('amor'):
|
||||
NICOS_CACHE_DIR = '/home/amor/nicosdata/amor/cache/'
|
||||
NICOS_CACHE_DIR = '/home/data/nicosdata/cache/'
|
||||
GREP = '/usr/bin/grep "value"'
|
||||
else:
|
||||
NICOS_CACHE_DIR = None
|
||||
@@ -74,25 +74,32 @@ class AmorHeader:
|
||||
logging.warning(f" parameter {key} not found, relpace by zero")
|
||||
return dtype(0)
|
||||
|
||||
def get_hdf_single_entry(self, path):
|
||||
if not np.shape(self.hdf['entry1/title']):
|
||||
return self.hdf[path][()].decode('utf-8')
|
||||
else:
|
||||
# format until 2025
|
||||
return self.hdf[path][0].decode('utf-8')
|
||||
|
||||
def read_header_info(self):
|
||||
# read general information and first data set
|
||||
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')
|
||||
title = self.get_hdf_single_entry('entry1/title')
|
||||
proposal_id = self.get_hdf_single_entry('entry1/proposal_id')
|
||||
user_name = self.get_hdf_single_entry('entry1/user/name')
|
||||
user_affiliation = 'unknown'
|
||||
user_email = self.hdf['entry1/user/email'][0].decode('utf-8')
|
||||
user_email = self.get_hdf_single_entry('entry1/user/email')
|
||||
user_orcid = None
|
||||
sampleName = self.hdf['entry1/sample/name'][0].decode('utf-8')
|
||||
model = self.hdf['entry1/sample/model'][0].decode('utf-8')
|
||||
if 'stack:' in model:
|
||||
sampleName = self.get_hdf_single_entry('entry1/sample/name')
|
||||
instrumentName = 'Amor'
|
||||
source = self.get_hdf_single_entry('entry1/Amor/source/name')
|
||||
sourceProbe = 'neutron'
|
||||
model = self.get_hdf_single_entry('entry1/sample/model')
|
||||
if 'stack' in model:
|
||||
import yaml
|
||||
model = yaml.safe_load(model)
|
||||
else:
|
||||
model = dict(stack=model)
|
||||
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_time = self.get_hdf_single_entry('entry1/start_time')
|
||||
# extract start time as unix time, adding UTC offset of 1h to time string
|
||||
start_date = datetime.fromisoformat(start_time)
|
||||
self.fileDate = start_date.replace(tzinfo=AMOR_LOCAL_TIMEZONE)
|
||||
@@ -113,11 +120,23 @@ class AmorHeader:
|
||||
facility=source,
|
||||
proposalID=proposal_id
|
||||
)
|
||||
if model['stack'] == '':
|
||||
om = None
|
||||
else:
|
||||
om = SampleModel.from_dict(model)
|
||||
self.sample = fileio.Sample(
|
||||
name=sampleName,
|
||||
model=SampleModel.from_dict(model),
|
||||
sample_parameters=None,
|
||||
model=om,
|
||||
sample_parameters={},
|
||||
)
|
||||
# while event times are not evaluated, use average_value reported in file for SEE
|
||||
if self.hdf['entry1/sample'].get('temperature', None) is not None \
|
||||
and float(self.hdf['entry1/sample/temperature/average_value'][0])>0:
|
||||
self.sample.sample_parameters['temperature'] = fileio.Value(
|
||||
float(self.hdf['entry1/sample/temperature/average_value'][0]), unit='K')
|
||||
if self.hdf['entry1/sample'].get('magnetic_field', None) is not None:
|
||||
self.sample.sample_parameters['magnetic_field'] = fileio.Value(
|
||||
float(self.hdf['entry1/sample/magnetic_field/average_value'][0]), unit='T')
|
||||
|
||||
def read_instrument_configuration(self):
|
||||
chopperSeparation = float(np.take(self.hdf['entry1/Amor/chopper/pair_separation'], 0))
|
||||
@@ -163,9 +182,11 @@ class AmorHeader:
|
||||
round(mu+kap+kad+0.5*div, 3),
|
||||
'deg'),
|
||||
wavelength = fileio.ValueRange(const.lamdaCut, const.lamdaMax, 'angstrom'),
|
||||
#polarization = fileio.Polarization.unpolarized,
|
||||
polarization = fileio.Polarization(polarizationConfig)
|
||||
)
|
||||
self.instrument_settings.qz = fileio.ValueRange(round(4*np.pi*np.sin(np.deg2rad(mu+kap+kad-0.5*div))/const.lamdaMax, 3),
|
||||
round(4*np.pi*np.sin(np.deg2rad(mu+kap+kad+0.5*div))/const.lamdaCut, 3),
|
||||
'1/angstrom')
|
||||
self.instrument_settings.mu = fileio.Value(
|
||||
round(mu, 3),
|
||||
'deg',
|
||||
|
||||
@@ -25,8 +25,8 @@ except ImportError:
|
||||
|
||||
|
||||
KAFKA_BROKER = 'linkafka01.psi.ch:9092'
|
||||
AMOR_EVENTS = 'amor_ev44'
|
||||
AMOR_NICOS = 'AMOR_nicosForwarder'
|
||||
AMOR_EVENTS = 'amor_detector'
|
||||
AMOR_NICOS = 'amor_nicosForwarder'
|
||||
|
||||
class KafkaFrozenData:
|
||||
"""
|
||||
|
||||
@@ -44,9 +44,9 @@ from .projection import TofZProjection, YZProjection
|
||||
|
||||
KAFKA_BROKER = 'linkafka01.psi.ch:9092'
|
||||
KAFKA_TOPICS = {
|
||||
'histogram': 'AMOR_histograms',
|
||||
'response': 'AMOR_histResponse',
|
||||
'command': 'AMOR_histCommands'
|
||||
'histogram': 'amor_histograms',
|
||||
'response': 'amor_histResponse',
|
||||
'command': 'amor_histCommands'
|
||||
}
|
||||
|
||||
def ktime():
|
||||
|
||||
@@ -6,6 +6,8 @@ import os
|
||||
import numpy as np
|
||||
from typing import List, Optional
|
||||
|
||||
from orsopy import fileio
|
||||
|
||||
from .event_data_types import EventDatasetProtocol
|
||||
from .header import Header
|
||||
from .options import NormalisationMethod
|
||||
@@ -24,7 +26,6 @@ class LZNormalisation:
|
||||
detZ_e = reference.data.events.detZ
|
||||
self.monitor = np.sum(reference.data.pulses.monitor)
|
||||
norm_lz, _, _ = np.histogram2d(lamda_e, detZ_e, bins=(grid.lamda(), grid.z()))
|
||||
norm_lz = np.where(norm_lz>2, norm_lz, np.nan)
|
||||
if normalisationMethod==NormalisationMethod.direct_beam:
|
||||
self.norm = np.flip(norm_lz, 1)
|
||||
else:
|
||||
@@ -47,7 +48,7 @@ class LZNormalisation:
|
||||
self = super().__new__(cls)
|
||||
with open(filename, 'rb') as fh:
|
||||
hash = str(np.load(fh, allow_pickle=True))
|
||||
self.file_list = np.load(fh, allow_pickle=True)
|
||||
self.file_list = np.load(fh, allow_pickle=True).tolist()
|
||||
self.angle = np.load(fh, allow_pickle=True)
|
||||
self.norm = np.load(fh, allow_pickle=True)
|
||||
self.monitor = np.load(fh, allow_pickle=True)
|
||||
@@ -99,4 +100,33 @@ class LZNormalisation:
|
||||
np.save(fh, self.monitor, allow_pickle=False)
|
||||
|
||||
def update_header(self, header:Header):
|
||||
header.measurement_additional_files = self.file_list
|
||||
header.measurement_additional_files = [fileio.File(file=os.path.basename(entry)) for entry in self.file_list]
|
||||
|
||||
def smooth(self, width):
|
||||
logging.debug(f'apply convolution with gaussian along lambda axis to smooth norm, sigma={width}')
|
||||
try:
|
||||
from scipy.signal import fftconvolve
|
||||
except ImportError:
|
||||
self._smooth_slow(width)
|
||||
kx = np.arange(self.norm.shape[0])-self.norm.shape[0]/2.
|
||||
kernel = np.exp(-0.5*kx**2/width**2)
|
||||
kernel/=kernel.sum()
|
||||
kernel = kernel[:, np.newaxis]*np.ones(self.norm.shape[1])[np.newaxis, :]
|
||||
unorm = np.where(np.isnan(self.norm), 0., self.norm)
|
||||
nnorm = fftconvolve(unorm, kernel, mode='same', axes=0)
|
||||
nnorm[np.isnan(self.norm)] = np.nan
|
||||
self.norm = nnorm
|
||||
|
||||
def _smooth_slow(self, width):
|
||||
# like smooth but using numpy buildin slow convolve
|
||||
nnorm = np.zeros_like(self.norm)
|
||||
|
||||
kx = np.arange(self.norm.shape[0])-self.norm.shape[0]/2.
|
||||
kernel = np.exp(-0.5*kx**2/width**2)
|
||||
kernel/=kernel.sum()
|
||||
unorm = np.where(np.isnan(self.norm), 0., self.norm)
|
||||
|
||||
for row in range(self.norm.shape[1]):
|
||||
nnorm[:, row] = np.convolve(unorm[:, row], kernel, mode='same')
|
||||
nnorm[np.isnan(self.norm)] = np.nan
|
||||
self.norm = nnorm
|
||||
|
||||
138
eos/options.py
138
eos/options.py
@@ -21,6 +21,11 @@ except ImportError:
|
||||
# python <3.10 use Enum instead
|
||||
from enum import Enum as StrEnum
|
||||
|
||||
class InCallString(StrEnum):
|
||||
auto='auto'
|
||||
always='always'
|
||||
never='never'
|
||||
|
||||
@dataclass
|
||||
class CommandlineParameterConfig:
|
||||
argument: str # default parameter for command line resutign ins "--argument"
|
||||
@@ -28,6 +33,7 @@ class CommandlineParameterConfig:
|
||||
short_form: Optional[str] = None
|
||||
group: str = 'misc'
|
||||
priority: int = 0
|
||||
in_call_string: InCallString = InCallString.auto
|
||||
|
||||
def __gt__(self, other):
|
||||
"""
|
||||
@@ -90,6 +96,7 @@ class ArgParsable:
|
||||
typ = field.type
|
||||
if get_origin(typ) is list:
|
||||
args['nargs'] = '+'
|
||||
args['action'] = 'extend'
|
||||
typ = get_args(typ)[0]
|
||||
if get_origin(typ) is tuple:
|
||||
# tuple of items are put together during evaluation
|
||||
@@ -117,6 +124,7 @@ class ArgParsable:
|
||||
group=field.metadata.get('group', 'misc'),
|
||||
short_form=field.metadata.get('short', None),
|
||||
priority=field.metadata.get('priority', 0),
|
||||
in_call_string=field.metadata.get('in_call_string', InCallString.auto),
|
||||
))
|
||||
return output
|
||||
|
||||
@@ -168,6 +176,34 @@ class ArgParsable:
|
||||
inpargs[field.name] = value
|
||||
return cls(**inpargs)
|
||||
|
||||
def get_call_parameters(self, abbrv=True):
|
||||
"""
|
||||
Return a list of command line arguments that reproduce this config, do not add default parameters.
|
||||
"""
|
||||
output = []
|
||||
for arg in sorted(self.get_commandline_parameters()):
|
||||
if ((arg.in_call_string==InCallString.auto and self.is_default(arg.argument)) or
|
||||
arg.in_call_string==InCallString.never):
|
||||
# skip default arguments or arguments defined to never appear in call string
|
||||
continue
|
||||
if arg.short_form and abbrv:
|
||||
item = '-' + arg.short_form + ' '
|
||||
else:
|
||||
item = '--' + arg.argument + ' '
|
||||
if arg.add_argument_args.get('type', None) in [str, float, int]:
|
||||
nargs = arg.add_argument_args.get('nargs', None)
|
||||
if nargs is None:
|
||||
item += str(getattr(self, arg.argument))
|
||||
elif nargs=='+':
|
||||
# remove the default parameters, only show added ones
|
||||
ignore = len(arg.add_argument_args.get('default', []))
|
||||
item += ' '.join([str(pi) for pi in getattr(self, arg.argument)[ignore:]])
|
||||
else:
|
||||
item += ' '.join([str(pi) for pi in getattr(self, arg.argument)])
|
||||
# boolean flags only reach this point if they are non-default
|
||||
output.append((arg, item))
|
||||
return output
|
||||
|
||||
# definition of command line arguments
|
||||
|
||||
@dataclass
|
||||
@@ -178,6 +214,7 @@ class ReaderConfig(ArgParsable):
|
||||
'short': 'Y',
|
||||
'group': 'input data',
|
||||
'help': 'year the measurement was performed',
|
||||
'in_call_string': InCallString.always,
|
||||
},
|
||||
)
|
||||
rawPath: List[str] = field(
|
||||
@@ -302,7 +339,7 @@ class ExperimentConfig(ArgParsable):
|
||||
},
|
||||
)
|
||||
muOffset: Optional[float] = field(
|
||||
default=0,
|
||||
default=None,
|
||||
metadata={
|
||||
'short': 'm',
|
||||
'group': 'sample',
|
||||
@@ -381,6 +418,21 @@ class ReflectivityReductionConfig(ArgParsable):
|
||||
'priority': 90,
|
||||
'group': 'input data',
|
||||
'help': 'normalisation method: [o]verillumination, [u]nderillumination, [d]irect_beam'})
|
||||
normalizationFilter: float = field(
|
||||
default=-1,
|
||||
metadata={
|
||||
'group': 'input data',
|
||||
'help': 'minimum normalization counts in lambda-theta bin to use, else filter'})
|
||||
normAngleFilter: float = field(
|
||||
default=-1,
|
||||
metadata={
|
||||
'group': 'input data',
|
||||
'help': 'minimum normalization counts total thetat bin to use, else filter'})
|
||||
normalizationSmoothing: float = field(
|
||||
default=0,
|
||||
metadata={
|
||||
'group': 'input data',
|
||||
'help': 'apply convolution on lambda axes to smooth the normalization data, useful for low statistics'})
|
||||
scale: List[float] = field(
|
||||
default_factory=lambda: [1.],
|
||||
metadata={
|
||||
@@ -404,7 +456,7 @@ class ReflectivityReductionConfig(ArgParsable):
|
||||
'group': 'input data',
|
||||
'help': 'File with R(q_z) curve to be subtracted (in .Rqz.ort format)'})
|
||||
normalisationFileIdentifier: Optional[List[str]] = field(
|
||||
default_factory=lambda: [None],
|
||||
default=None,
|
||||
metadata={
|
||||
'short': 'n',
|
||||
'priority': 90,
|
||||
@@ -526,85 +578,19 @@ class ReflectivityConfig:
|
||||
|
||||
def call_string(self):
|
||||
base = 'eos'
|
||||
|
||||
inpt = ''
|
||||
if self.reader.year:
|
||||
inpt += f' -Y {self.reader.year}'
|
||||
else:
|
||||
inpt += f' -Y {datetime.now().year}'
|
||||
if np.shape(self.reader.rawPath)[0] == 1:
|
||||
inpt += f' --rawPath {self.reader.rawPath}'
|
||||
if self.reduction.subtract:
|
||||
inpt += f' -subtract {self.reduction.subtract}'
|
||||
if self.reduction.normalisationFileIdentifier:
|
||||
inpt += f' -n {" ".join(self.reduction.normalisationFileIdentifier)}'
|
||||
if self.reduction.fileIdentifier:
|
||||
inpt += f' -f {" ".join(self.reduction.fileIdentifier)}'
|
||||
|
||||
otpt = ''
|
||||
if self.reduction.qResolution:
|
||||
otpt += f' -r {self.reduction.qResolution}'
|
||||
if self.output.outputPath != '.':
|
||||
inpt += f' --outputdPath {self.output.outputPath}'
|
||||
if self.output.outputName:
|
||||
otpt += f' -o {self.output.outputName}'
|
||||
if self.output.outputFormats != ['Rqz.ort']:
|
||||
otpt += f' -of {" ".join(self.output.outputFormats)}'
|
||||
|
||||
mask = ''
|
||||
call_parameters = self.reader.get_call_parameters()
|
||||
call_parameters += self.output.get_call_parameters()
|
||||
call_parameters += self.reduction.get_call_parameters()
|
||||
call_parameters += self.experiment.get_call_parameters()
|
||||
|
||||
mask += f' -y {" ".join(str(ii) for ii in self.experiment.yRange)}'
|
||||
mask += f' -l {" ".join(str(ff) for ff in self.experiment.lambdaRange)}'
|
||||
mask += f' -t {" ".join(str(ff) for ff in self.reduction.thetaRange)}'
|
||||
mask += f' -T {" ".join(str(ff) for ff in self.reduction.thetaRangeR)}'
|
||||
mask += f' -q {" ".join(str(ff) for ff in self.reduction.qzRange)}'
|
||||
call_parameters.sort()
|
||||
|
||||
para = ''
|
||||
# TODO: Check if we want these parameters for defaults
|
||||
para += f' --chopperPhase {self.experiment.chopperPhase}'
|
||||
para += f' --chopperPhaseOffset {self.experiment.chopperPhaseOffset}'
|
||||
if self.experiment.mu:
|
||||
para += f' --mu {self.experiment.mu}'
|
||||
elif self.experiment.muOffset:
|
||||
para += f' --muOffset {self.experiment.muOffset}'
|
||||
if self.experiment.nu:
|
||||
para += f' --nu {self.experiment.nu}'
|
||||
cpout = f'{base} ' + ' '.join([cp[1] for cp in call_parameters])
|
||||
|
||||
modl = ''
|
||||
if self.experiment.sampleModel:
|
||||
modl += f" --sampleModel '{self.experiment.sampleModel}'"
|
||||
|
||||
acts = ''
|
||||
if self.reduction.autoscale:
|
||||
acts += f' --autoscale {" ".join(str(ff) for ff in self.reduction.autoscale)}'
|
||||
# TODO: Check if should be shown if not default
|
||||
acts += f' --scale {self.reduction.scale}'
|
||||
if self.reduction.timeSlize:
|
||||
acts += f' --timeSlize {" ".join(str(ff) for ff in self.reduction.timeSlize)}'
|
||||
|
||||
mlst = base + inpt + otpt
|
||||
if mask:
|
||||
mlst += mask
|
||||
if para:
|
||||
mlst += para
|
||||
if acts:
|
||||
mlst += acts
|
||||
if modl:
|
||||
mlst += modl
|
||||
|
||||
if len(mlst) > 70:
|
||||
mlst = base + ' ' + inpt + ' ' + otpt
|
||||
if mask:
|
||||
mlst += ' ' + mask
|
||||
if para:
|
||||
mlst += ' ' + para
|
||||
if acts:
|
||||
mlst += ' ' + acts
|
||||
if modl:
|
||||
mlst += ' ' + modl
|
||||
|
||||
logging.debug(f'Argument list build in EOSConfig.call_string: {mlst}')
|
||||
return mlst
|
||||
logging.debug(f'Argument list build in EOSConfig.call_string: {cpout}')
|
||||
return cpout
|
||||
|
||||
class E2HPlotSelection(StrEnum):
|
||||
All = 'all'
|
||||
|
||||
@@ -85,12 +85,19 @@ class ProjectedReflectivity:
|
||||
self.R -= R
|
||||
self.dR = np.sqrt(self.dR**2+dR**2)
|
||||
|
||||
def plot(self, **kwargs):
|
||||
from matplotlib import pyplot as plt
|
||||
plt.errorbar(self.Q, self.R, yerr=self.dR, **kwargs)
|
||||
plt.yscale('log')
|
||||
plt.xlabel('Q / Å$^{-1}$')
|
||||
plt.ylabel('R')
|
||||
|
||||
class LZProjection(ProjectionInterface):
|
||||
grid: LZGrid
|
||||
lamda: np.ndarray
|
||||
alphaF: np.ndarray
|
||||
is_normalized: bool
|
||||
angle: float
|
||||
|
||||
data: np.recarray
|
||||
_dtype = np.dtype([
|
||||
@@ -107,6 +114,7 @@ class LZProjection(ProjectionInterface):
|
||||
def __init__(self, tthh: float, grid: LZGrid):
|
||||
self.grid = grid
|
||||
self.is_normalized = False
|
||||
self.angle = tthh
|
||||
|
||||
alphaF_z = tthh + Detector.delta_z
|
||||
lamda_l = self.grid.lamda()
|
||||
@@ -168,9 +176,12 @@ class LZProjection(ProjectionInterface):
|
||||
self.data.mask &= self.lamda>=lamda_range[0]
|
||||
self.data.mask &= self.lamda<=lamda_range[1]
|
||||
|
||||
def apply_norm_mask(self, norm: LZNormalisation):
|
||||
def apply_norm_mask(self, norm: LZNormalisation, min_norm=-1, min_theta=-1):
|
||||
# Mask points where normliazation is nan
|
||||
self.data.mask &= np.logical_not(np.isnan(norm.norm))
|
||||
self.data.mask &= np.logical_not(np.isnan(norm.norm))&(norm.norm>min_norm)
|
||||
if min_theta>0:
|
||||
thsum = np.nansum(norm.norm, axis=0)
|
||||
self.data.mask &= (thsum>min_theta)[np.newaxis, :]
|
||||
|
||||
def project(self, dataset: EventDatasetProtocol, monitor: float):
|
||||
"""
|
||||
@@ -190,6 +201,7 @@ class LZProjection(ProjectionInterface):
|
||||
self.data[:] = 0
|
||||
self.data.mask = True
|
||||
self.monitor = 0.
|
||||
self.norm_monitor = 1.
|
||||
|
||||
@property
|
||||
def I(self):
|
||||
@@ -199,22 +211,25 @@ class LZProjection(ProjectionInterface):
|
||||
|
||||
def calc_error(self):
|
||||
# calculate error bars for resulting intensity after normalization
|
||||
self.data.err = self.data.ref * np.sqrt( 1./(self.data.I+.1) + 1./self.data.norm )
|
||||
self.data.err = self.data.ref * np.sqrt( 1./(self.data.I+.1) + 1./(self.data.norm+0.1) )
|
||||
|
||||
def normalize_over_illuminated(self, norm: LZNormalisation):
|
||||
"""
|
||||
Normalize the dataaset and take into account a difference in
|
||||
detector angle for measurement and reference.
|
||||
"""
|
||||
logging.debug(f' correcting for incident angle difference from norm {norm.angle} to data {self.angle}')
|
||||
norm_lz = norm.norm
|
||||
thetaN_z = Detector.delta_z+norm.angle
|
||||
thetaN_lz = np.ones_like(norm_lz)*thetaN_z
|
||||
thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan)
|
||||
self.data.mask &= (np.absolute(thetaN_lz)>5e-3)
|
||||
ref_lz = (self.data.I*np.absolute(thetaN_lz))/(norm_lz*np.absolute(self.alphaF))
|
||||
delta_lz = np.ones_like(norm_lz)*Detector.delta_z
|
||||
# do not perform gravity correction for footprint, would require norm detector distance that is unknown here
|
||||
fp_corr_lz = np.where(np.absolute(delta_lz+norm.angle)>5e-3,
|
||||
(delta_lz+self.angle)/(delta_lz+norm.angle), np.nan)
|
||||
self.data.mask &= np.logical_not(np.isnan(fp_corr_lz))
|
||||
self.data.norm = norm_lz*fp_corr_lz
|
||||
self.norm_monitor = norm.monitor
|
||||
ref_lz = self.data.I/np.where(self.data.norm>0, self.data.norm, np.nan)
|
||||
ref_lz *= norm.monitor/self.monitor
|
||||
ref_lz[np.logical_not(self.data.mask)] = np.nan
|
||||
self.data.norm = norm_lz
|
||||
self.data.ref = ref_lz
|
||||
self.calc_error()
|
||||
self.is_normalized = True
|
||||
@@ -234,6 +249,7 @@ class LZProjection(ProjectionInterface):
|
||||
raise ValueError("Dataset needs to be normalized, first")
|
||||
self.data.ref *= factor
|
||||
self.data.err *= factor
|
||||
self.norm_monitor /= factor
|
||||
|
||||
def project_on_qz(self):
|
||||
if not self.is_normalized:
|
||||
@@ -241,29 +257,25 @@ class LZProjection(ProjectionInterface):
|
||||
q_q = self.grid.q()
|
||||
weights_lzf = self.data.norm[self.data.mask]
|
||||
q_lzf = self.data.qz[self.data.mask]
|
||||
R_lzf = self.data.ref[self.data.mask]
|
||||
dR_lzf = self.data.err[self.data.mask]
|
||||
I_lzf = self.data.I[self.data.mask]
|
||||
dq_lzf = self.data.res[self.data.mask]
|
||||
|
||||
N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf )[0]
|
||||
# get number of grid points contributing to a bin, filter points with no contribution
|
||||
N_q = np.histogram(q_lzf, bins = q_q)[0]
|
||||
N_q = np.where(N_q > 0, N_q, np.nan)
|
||||
fltr = N_q>0
|
||||
|
||||
R_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf * R_lzf )[0]
|
||||
R_q = R_q / N_q
|
||||
|
||||
dR_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dR_lzf)**2 )[0]
|
||||
dR_q = np.sqrt( dR_q ) / N_q
|
||||
|
||||
# TODO: different error propagations for dR and dq!
|
||||
# this is what should work:
|
||||
#dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0]
|
||||
#dq_q = np.sqrt( dq_q ) / N_q
|
||||
# and this actually works:
|
||||
N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf**2 )[0]
|
||||
N_q = np.where(N_q > 0, N_q, np.nan)
|
||||
dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0]
|
||||
dq_q = np.sqrt( dq_q / N_q )
|
||||
|
||||
# calculate sum of all normalization weights per bin
|
||||
W_q = np.maximum(np.histogram(q_lzf, bins = q_q, weights = weights_lzf)[0], 1e-10)
|
||||
# calculate sum of all dataset counts per bin
|
||||
I_q = np.histogram(q_lzf, bins = q_q, weights = I_lzf)[0]
|
||||
# normlaize dataaset by normalization counts and scale by monitor
|
||||
R_q = np.where(fltr, I_q*self.norm_monitor/self.monitor / W_q, np.nan)
|
||||
# error as squar-root of counts and sqrt from normalization (dR/R = sqrt( (dI/I)² + (dW/W)²)
|
||||
dR_q = np.where(fltr, R_q*(np.sqrt(1./(I_q+0.1)+ 1./(W_q+0.1))), np.nan)
|
||||
# q-resolution is the weighted sum of individual points q-resolutions
|
||||
dq_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf * dq_lzf )[0]
|
||||
dq_q = np.where(fltr, dq_q/W_q, np.nan)
|
||||
return ProjectedReflectivity(R_q, dR_q, (q_q[1:]+q_q[:-1])/2., dq_q)
|
||||
|
||||
def plot(self, **kwargs):
|
||||
|
||||
@@ -236,7 +236,7 @@ class E2HReduction:
|
||||
self.event_actions(self.dataset)
|
||||
self.dataset.update_header(self.header)
|
||||
|
||||
self.header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1],
|
||||
self.header.measurement_data_files.append(fileio.File(file=os.path.basename(fileName),
|
||||
timestamp=self.dataset.fileDate))
|
||||
|
||||
def add_data(self):
|
||||
|
||||
@@ -85,6 +85,9 @@ class ReflectivityReduction:
|
||||
else:
|
||||
self.norm = LZNormalisation.unity(self.grid)
|
||||
|
||||
if self.config.reduction.normalizationSmoothing:
|
||||
self.norm.smooth(self.config.reduction.normalizationSmoothing)
|
||||
|
||||
# load R(q_z) curve to be subtracted:
|
||||
if self.config.reduction.subtract:
|
||||
self.sq_q, self.sR_q, self.sdR_q, self.sFileName = self.loadRqz(self.config.reduction.subtract)
|
||||
@@ -101,6 +104,16 @@ class ReflectivityReduction:
|
||||
self.read_file_block(i, short_notation)
|
||||
|
||||
# output
|
||||
if self.config.output.is_default('outputName'):
|
||||
import datetime
|
||||
_date = datetime.datetime.now().replace(microsecond=0).isoformat()
|
||||
if self.header.sample.name:
|
||||
_sampleName = self.header.sample.name.replace(' ', '_')
|
||||
else:
|
||||
_sampleName = 'unknown'
|
||||
_mu = int(self.dataset.geometry.mu * 3)
|
||||
self.config.output.outputName = f'{_sampleName}_{_mu:03}_{_date}'
|
||||
|
||||
logging.warning('output:')
|
||||
|
||||
if 'Rqz.ort' in self.config.output.outputFormats:
|
||||
@@ -146,7 +159,7 @@ class ReflectivityReduction:
|
||||
self.dataset.append(di)
|
||||
|
||||
for fileName in file_list:
|
||||
self.header.measurement_data_files.append(fileio.File( file=fileName.split('/')[-1],
|
||||
self.header.measurement_data_files.append(fileio.File( file=os.path.basename(fileName),
|
||||
timestamp=self.dataset.fileDate))
|
||||
|
||||
|
||||
@@ -364,7 +377,7 @@ class ReflectivityReduction:
|
||||
if reference.data.events.shape[0] > 1e6:
|
||||
self.norm.safe(n_path, self.normevent_actions.action_hash())
|
||||
|
||||
self.header.measurement_additional_files = self.norm.file_list
|
||||
self.norm.update_header(self.header)
|
||||
self.header.reduction.corrections.append('normalisation with \'additional files\'')
|
||||
|
||||
def project_on_lz(self, dataset=None):
|
||||
@@ -390,7 +403,8 @@ class ReflectivityReduction:
|
||||
|
||||
proj.apply_lamda_mask(self.config.experiment.lambdaRange)
|
||||
|
||||
proj.apply_norm_mask(self.norm)
|
||||
proj.apply_norm_mask(self.norm, min_norm=self.config.reduction.normalizationFilter,
|
||||
min_theta=self.config.reduction.normAngleFilter)
|
||||
|
||||
proj.project(dataset, self.monitor)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user