Merge pull request #9 from jochenstahn/new_filewriter

New filewriter

eos/__init__.py

@@ -2,5 +2,5 @@
 Package to handle data redction at AMOR instrument to be used by __main__.py script.
 """
 
-__version__ = '3.0.3'
-__date__    = '2025-10-28'
+__version__ = '3.1.0'
+__date__    = '2026-02-20'

eos/event_analysis.py

@@ -25,8 +25,15 @@ class ExtractWalltime(EventDataAction):
         dataset.data.events = new_events
 
 class MergeFrames(EventDataAction):
+    def __init__(self, lamdaCut=None):
+        self.lamdaCut=lamdaCut
+
     def perform_action(self, dataset: EventDatasetProtocol)->None:
-        tofCut = const.lamdaCut*dataset.geometry.chopperDetectorDistance/const.hdm*1e-13
+        if self.lamdaCut is None:
+            lamdaCut = const.lamdaCut
+        else:
+            lamdaCut = self.lamdaCut
+        tofCut = lamdaCut*dataset.geometry.chopperDetectorDistance/const.hdm*1e-13
         total_offset = (tofCut +
                         dataset.timing.tau * (dataset.timing.ch1TriggerPhase + dataset.timing.chopperPhase/2)/180)
         dataset.data.events.tof = merge_frames(dataset.data.events.tof, tofCut, dataset.timing.tau, total_offset)

eos/file_reader.py

@@ -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',
@@ -223,7 +244,6 @@ class AmorEventData(AmorHeader):
 
     def __init__(self, fileName:Union[str, h5py.File, BinaryIO], first_index:int=0, max_events:int=100_000_000):
         if type(fileName) is str:
-            #logging.warning(f'    opening file {fileName}')
             logging.warning(f'    {fileName.split("/")[-1]}')
             self.file_list = [fileName]
             hdf = h5py.File(fileName, 'r', swmr=True)

eos/instrument.py

@@ -65,9 +65,14 @@ class LZGrid:
     def qzRange(self):
         return self._qzRange
 
-    def __init__(self, qResolution, qzRange):
+    def __init__(self, qResolution, qzRange, lambda_overwrite=None):
         self._qResolution = qResolution
         self._qzRange = qzRange
+        if lambda_overwrite is None:
+            self.lamdaMax = const.lamdaMax
+            self.lamdaCut = const.lamdaCut
+        else:
+            self.lamdaCut, self.lamdaMax = lambda_overwrite
 
     @property
     @cache
@@ -92,8 +97,8 @@ class LZGrid:
 
     @cache
     def lamda(self):
-        lamdaMax = 16
-        lamdaMin = const.lamdaCut
+        lamdaMax = self.lamdaMax
+        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
 

eos/kafka_events.py

@@ -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:
     """

eos/kafka_serializer.py

@@ -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():

eos/normalization.py

@@ -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

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'

eos/projection.py

@@ -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):
@@ -298,7 +310,7 @@ class LZProjection(ProjectionInterface):
                 plt.colorbar(label='I / cpm')
         plt.xlabel('$\\lambda$ / $\\AA$')
         plt.ylabel('$\\Theta$ / °')
-        plt.xlim(3., 12.)
+        plt.xlim(self.lamda[0,0], self.lamda[-1,0])
         af = self.alphaF[self.data.mask]
         plt.ylim(af.min(), af.max())
         plt.title('Wavelength vs. Reflection Angle')

eos/reduction_e2h.py

@@ -51,6 +51,7 @@ class E2HReduction:
         self.plot_kwds = {}
         plt.rcParams.update({'font.size': self.config.reduction.fontsize})
 
+        self.overwrite = eh.ApplyParameterOverwrites(self.config.experiment) # some actions use instrument parameters, change before that
         if self.config.reduction.update:
             # live update implies plotting
             self.config.reduction.show_plot = True
@@ -58,6 +59,9 @@ class E2HReduction:
         if self.config.reduction.plot==E2HPlotSelection.Raw:
             # Raw implies fast caculations
             self.config.reduction.fast = True
+        if not self.config.experiment.is_default('lambdaRange'):
+            # filtering wavelength requires frame analysis
+            self.config.reduction.fast = False
 
         if not self.config.reduction.fast or self.config.reduction.plot in NEEDS_LAMDA:
             from . import event_analysis as ea
@@ -65,6 +69,7 @@ class E2HReduction:
         # Actions on datasets not used for normalization
         self.event_actions = eh.ApplyPhaseOffset(self.config.experiment.chopperPhaseOffset)
         if not self.config.reduction.fast:
+            self.event_actions |= self.overwrite
             self.event_actions |= eh.CorrectChopperPhase()
             self.event_actions |= ea.ExtractWalltime()
         else:
@@ -84,8 +89,8 @@ class E2HReduction:
                 # perform corrections for tof if not fast mode
                 self.event_actions |= eh.TofTimeCorrection(self.config.experiment.incidentAngle==IncidentAngle.alphaF)
         # select needed actions in depenence of plots
-        if self.config.reduction.plot in NEEDS_LAMDA:
-            self.event_actions |= ea.MergeFrames()
+        if self.config.reduction.plot in NEEDS_LAMDA or not self.config.experiment.is_default('lambdaRange'):
+            self.event_actions |= ea.MergeFrames(lamdaCut=self.config.experiment.lambdaRange[0])
             self.event_actions |= ea.AnalyzePixelIDs(self.config.experiment.yRange)
             self.event_actions |= eh.TofTimeCorrection(self.config.experiment.incidentAngle==IncidentAngle.alphaF)
             self.event_actions |= ea.CalculateWavelength(self.config.experiment.lambdaRange)
@@ -93,8 +98,8 @@ class E2HReduction:
 
         # plot dependant options
         if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]:
-            self.grid = LZGrid(0.05, [0.0, 0.25])
-            self.grid.dldl = 0.05
+            self.grid = LZGrid(0.05, [0.0, 0.25], lambda_overwrite=self.config.experiment.lambdaRange)
+            self.grid.dldl = 0.01
 
         if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.Raw,
                                           E2HPlotSelection.LT, E2HPlotSelection.YT,
@@ -148,6 +153,7 @@ class E2HReduction:
 
     def prepare_graphs(self):
         last_file_header = AmorHeader(self.file_list[-1])
+        self.overwrite.perform_action(last_file_header)
         tthh  = last_file_header.geometry.nu - last_file_header.geometry.mu
 
         if not self.config.reduction.is_default('thetaRangeR'):
@@ -230,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):

eos/reduction_reflectivity.py

@@ -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:
@@ -117,7 +130,7 @@ class ReflectivityReduction:
             plt.show()
 
     def read_file_block(self, i, short_notation):
-        logging.warning('reading input:')
+        logging.warning('input:')
         file_list = self.path_resolver.resolve(short_notation)
 
         self.header.measurement_data_files = []
@@ -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))
 
 
@@ -159,7 +172,7 @@ class ReflectivityReduction:
 
     def analyze_unsliced(self, i):
         self.monitor = self.dataset.data.pulses.monitor.sum()
-        logging.warning(f'    monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.config.experiment.monitorType]}')
+        logging.info(f'    monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.config.experiment.monitorType]}')
 
         proj:LZProjection = self.project_on_lz()
         try:
@@ -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)
 

nicos_config.md

@@ -37,6 +37,7 @@ hist_tofz = device('nicos_sinq.devices.just_bin_it.JustBinItImage',
 ```
 
 These images have then to be set in the detector configuration as _images_ items:
+
 ```                      
 images=['hist_yz', 'hist_tofz'],
-```
+```