Add unit tests for filterinc capability and automatic spin state splitting

Co-authored-by: bruhn_b <basil.bruhn@psi.ch>
Reviewed-on: #1
Co-authored-by: Artur Glavic <artur.glavic@psi.ch>
Co-committed-by: Artur Glavic <artur.glavic@psi.ch>

.github/workflows/release.yml

@@ -22,7 +22,35 @@ on:
           - all_incl_release
 
 jobs:
+  test:
+
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ['3.8', '3.9', '3.10', '3.12']
+      fail-fast: false
+
+    steps:
+      - name: Checkout LFS objects
+        run: git clone https://${{secrets.GITHUB_TOKEN}}@gitea.psi.ch/${{ github.repository }}.git .
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install pytest
+          pip install -r requirements.txt
+
+      - name: Test with pytest
+        run: |
+          python -m pytest tests
+  
+
   build-ubuntu-latest:
+    needs: [test]
     runs-on: ubuntu-latest
     if: ${{ (github.event_name != 'workflow_dispatch') || (contains(fromJson('["all", "linux", "all_incl_release"]'), github.event.inputs.build-items)) }}
     permissions:
@@ -31,30 +59,28 @@ jobs:
     - uses: actions/checkout@v4
     - uses: actions/setup-python@v5
       with:
-        python-version: '3.11'
+        python-version: '3.12'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
         pip install build
         pip install -r requirements.txt
+        pip install wheel build twine
     - name: Build PyPI package
       run: |
         python3 -m build
-    - name: Archive distribution
-      uses: actions/upload-artifact@v4
-      with:
-        name: linux-dist
-        path: |
-          dist/*.tar.gz
+#    - name: Archive distribution
+#      uses: actions/upload-artifact@v3
+#      with:
+#        name: linux-dist
+#        path: |
+#          dist/*.tar.gz
     - name: Upload to PyPI
-      #if: github.event_name != 'workflow_dispatch'
-      uses: pypa/gh-action-pypi-publish@release/v1
-      with:
-      #  user: __token__
-      #  password: ${{ secrets.PYPI_TOKEN }}
-        skip-existing: true
+      run: |
+        twine upload dist/* -u __token__ -p ${{ secrets.PYPI_TOKEN }} --skip-existing
 
   build-windows:
+    needs: [test]
     runs-on: windows-latest
     if: ${{ (github.event_name != 'workflow_dispatch') || (contains(fromJson('["all", "windows", "all_incl_release"]'), github.event.inputs.build-items)) }}
 
@@ -74,7 +100,7 @@ jobs:
         cd dist\eos
         Compress-Archive -Path .\* -Destination ..\..\eos.zip
     - name: Archive distribution
-      uses: actions/upload-artifact@v4
+      uses: actions/upload-artifact@v3
       with:
         name: windows-dist
         path: |
@@ -89,10 +115,10 @@ jobs:
       with:
         fetch-depth: 0
         fetch-tags: true
-    - uses: actions/download-artifact@v4
+    - uses: actions/download-artifact@v3
       with:
         name: linux-dist
-    - uses: actions/download-artifact@v4
+    - uses: actions/download-artifact@v3
       with:
         name: windows-dist
     - name: get latest version tag

.github/workflows/unit_tests.yml

@@ -10,18 +10,18 @@ on:
   workflow_dispatch:
 
 jobs:
-  build:
+  test:
 
-    runs-on: ubuntu-22.04
+    runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ['3.8', '3.9', '3.10', '3.11', '3.12']
+        python-version: ['3.8', '3.9', '3.10', '3.12']
       fail-fast: false
 
     steps:
-      - uses: actions/checkout@v4
-        with:
-            lfs: 'true'
+      - name: Checkout LFS objects
+        run: git clone https://${{secrets.GITHUB_TOKEN}}@gitea.psi.ch/${{ github.repository }}.git .
+
       - name: Set up Python ${{ matrix.python-version }}
         uses: actions/setup-python@v5
         with:

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.0'
-__date__    = '2025-10-06'
+__version__ = '3.2.2'
+__date__    = '2026-02-27'

eos/__main__.py

@@ -8,30 +8,31 @@ Author: Jochen Stahn (algorithms, python draft),
 import logging
 
 # need to do absolute import here as pyinstaller requires it
-from eos.options import EOSConfig, ReaderConfig, ExperimentConfig, ReductionConfig, OutputConfig
+from eos.options import ReflectivityConfig, ReaderConfig, ExperimentConfig, ReflectivityReductionConfig, ReflectivityOutputConfig
 from eos.command_line import commandLineArgs
 from eos.logconfig import setup_logging, update_loglevel
 
+
 def main():
     setup_logging()
 
     # read command line arguments and generate classes holding configuration parameters
-    clas = commandLineArgs([ReaderConfig, ExperimentConfig, ReductionConfig, OutputConfig],
+    clas = commandLineArgs([ReaderConfig, ExperimentConfig, ReflectivityReductionConfig, ReflectivityOutputConfig],
                            'eos')
     update_loglevel(clas.verbose)
 
     reader_config = ReaderConfig.from_args(clas)
     experiment_config = ExperimentConfig.from_args(clas)
-    reduction_config = ReductionConfig.from_args(clas)
-    output_config = OutputConfig.from_args(clas)
-    config = EOSConfig(reader_config, experiment_config, reduction_config, output_config)
+    reduction_config = ReflectivityReductionConfig.from_args(clas)
+    output_config = ReflectivityOutputConfig.from_args(clas)
+    config = ReflectivityConfig(reader_config, experiment_config, reduction_config, output_config)
 
     logging.warning('######## eos - data reduction for Amor ########')
 
     # only import heavy module if sufficient command line parameters were provided
-    from eos.reduction import AmorReduction
+    from eos.reduction_reflectivity import ReflectivityReduction
     # Create reducer with these arguments
-    reducer = AmorReduction(config)
+    reducer = ReflectivityReduction(config)
     # Perform actual reduction
     reducer.reduce()
 

eos/command_line.py

@@ -1,10 +1,10 @@
 import argparse
 
-from typing import List
+from typing import List, Type
 from .options import ArgParsable
 
 
-def commandLineArgs(config_items: List[ArgParsable], program_name=None):
+def commandLineArgs(config_items: List[Type[ArgParsable]], program_name=None, extra_args=[]):
     """
     Process command line argument.
     The type of the default values is used for conversion and validation.
@@ -36,4 +36,7 @@ def commandLineArgs(config_items: List[ArgParsable], program_name=None):
                     f'--{cpc.argument}', **cpc.add_argument_args
                     )
 
+    for ma in extra_args:
+        clas.add_argument(**ma)
+
     return clas.parse_args()

eos/const.py

@@ -1,7 +1,11 @@
-"""
-Constants used in data reduction.
-"""
-
-hdm = 6.626176e-34/1.674928e-27  # h / m
-lamdaCut = 2.5  # Aa
-lamdaMax = 15.0 # Aa
+"""
+Constants used in data reduction.
+"""
+
+hdm = 6.626176e-34/1.674928e-27  # h / m
+lamdaCut = 2.5  # Aa
+lamdaMax = 15.0 # Aa
+
+polarizationConfigs = ['unpolarized', 'unpolarized', 'po', 'mo', 'op', 'pp', 'mp', 'om', 'pm', 'mm']
+polarizationLabels = ['undetermined', 'unpolarized', 'spin-up', 'spin-down', 'op',
+                      'up-up', 'down-up', 'om', 'up-down', 'down-down']

eos/e2h.py

@@ -0,0 +1,42 @@
+"""
+events2histogram vizualising data from Amor@SINQ, PSI
+
+Author: Jochen Stahn (algorithms, python draft),
+        Artur Glavic (structuring and optimisation of code)
+"""
+import logging
+
+# need to do absolute import here as pyinstaller requires it
+from eos.options import E2HConfig, ReaderConfig, ExperimentConfig, E2HReductionConfig
+from eos.command_line import commandLineArgs
+from eos.logconfig import setup_logging, update_loglevel
+
+
+def main():
+    setup_logging()
+    logging.getLogger('matplotlib').setLevel(logging.WARNING)
+
+    # read command line arguments and generate classes holding configuration parameters
+    clas = commandLineArgs([ReaderConfig, ExperimentConfig, E2HReductionConfig],
+                           'events2histogram')
+    update_loglevel(clas.verbose)
+
+    reader_config = ReaderConfig.from_args(clas)
+    experiment_config = ExperimentConfig.from_args(clas)
+    reduction_config = E2HReductionConfig.from_args(clas)
+    config = E2HConfig(reader_config, experiment_config, reduction_config)
+
+    logging.warning('######## events2histogram - data vizualization for Amor ########')
+    from eos.reduction_e2h import E2HReduction
+
+    # only import heavy module if sufficient command line parameters were provided
+    from eos.reduction_reflectivity import ReflectivityReduction
+    # Create reducer with these arguments
+    reducer = E2HReduction(config)
+    # Perform actual reduction
+    reducer.reduce()
+
+    logging.info('######## events2histogram - finished ########')
+
+if __name__ == '__main__':
+    main()

eos/event_analysis.py

@@ -9,7 +9,7 @@ from typing import Tuple
 
 from . import const
 from .event_data_types import EventDataAction, EventDatasetProtocol, append_fields, EVENT_BITMASKS
-from .helpers import filter_project_x, merge_frames, extract_walltime
+from .helpers import filter_project_x, merge_frames, extract_walltime, add_log_to_pulses
 from .instrument import Detector
 from .options import IncidentAngle
 from .header import Header
@@ -18,15 +18,22 @@ class ExtractWalltime(EventDataAction):
     def perform_action(self, dataset: EventDatasetProtocol) ->None:
         wallTime = extract_walltime(dataset.data.events.tof,
                                     dataset.data.packets.start_index,
-                                    dataset.data.packets.Time)
+                                    dataset.data.packets.time)
         logging.debug(f'        expending event stream by wallTime')
         new_events = append_fields(dataset.data.events, [('wallTime', wallTime.dtype)])
         new_events.wallTime = wallTime
         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)
@@ -117,5 +124,40 @@ class FilterQzRange(EventDataAction):
         if not 'qz' in dataset.data.events.dtype.names:
             raise ValueError("FilterQzRange requires dataset with qz values per events, perform WavelengthAndQ first")
 
-        if self.qzRange[1]<0.5:
-            d.events.mask += EVENT_BITMASKS["qRange"]*((self.qzRange[0]>d.events.qz) | (d.events.qz>self.qzRange[1]))
+        d.events.mask += EVENT_BITMASKS["qRange"]*((self.qzRange[0]>d.events.qz) | (d.events.qz>self.qzRange[1]))
+
+class FilterByLog(EventDataAction):
+
+    def __init__(self, filter_string, remove_switchpulse=False):
+        if filter_string.startswith('!'):
+            filter_string = filter_string[1:]
+            remove_switchpulse = True
+        self.filter_string = filter_string
+        self.remove_switchpulse = remove_switchpulse
+
+    def perform_action(self, dataset: EventDatasetProtocol) -> None:
+        filter_variable = None
+        # go through existing keys in reverse order of their length to insure a name containing another is used
+        existing_keys = list(dataset.data.device_logs.keys())
+        existing_keys.sort(key=lambda x: -len(x))
+        for key in existing_keys:
+            if key in self.filter_string:
+                filter_variable = key
+                break
+        if filter_variable is None:
+            logging.warning(f'          could not find suitable parameter to filter on in {self.filter_string}, '
+                            f'available parameters are: {list(sorted(existing_keys))}')
+            return
+        logging.debug(f'          using parameter {filter_variable} to apply filter {self.filter_string}')
+        if not filter_variable in EVENT_BITMASKS:
+            EVENT_BITMASKS[filter_variable] = max(EVENT_BITMASKS.values())*2
+        if not filter_variable in dataset.data.pulses.dtype.names:
+            # interpolate the parameter values for all existing pulses
+            add_log_to_pulses(filter_variable, dataset)
+        fltr_pulses = eval(self.filter_string, {filter_variable: dataset.data.pulses[filter_variable]})
+        if self.remove_switchpulse:
+            switched = fltr_pulses[:-1] & ~fltr_pulses[1:]
+            fltr_pulses[:-1] &= ~switched
+        goodTimeS = dataset.data.pulses.time[fltr_pulses]
+        filter_e = np.logical_not(np.isin(dataset.data.events.wallTime, goodTimeS))
+        dataset.data.events.mask += EVENT_BITMASKS[filter_variable]*filter_e

eos/event_data_types.py

@@ -1,8 +1,8 @@
 """
 Specify the data type and protocol used for event datasets.
 """
-from typing import List, Optional, Protocol, Tuple
-from dataclasses import dataclass
+from typing import Dict, List, Optional, Protocol, Tuple
+from dataclasses import dataclass, field
 from .header import Header
 from abc import ABC, abstractmethod
 from hashlib import sha256
@@ -31,9 +31,10 @@ class AmorTiming:
 
 # Structured datatypes used for event streams
 EVENT_TYPE = np.dtype([('tof', np.float64), ('pixelID', np.uint32), ('mask', np.int32)])
-PACKET_TYPE = np.dtype([('start_index', np.uint32), ('Time', np.int64)])
+PACKET_TYPE = np.dtype([('start_index', np.uint32), ('time', np.int64)])
 PULSE_TYPE = np.dtype([('time', np.int64), ('monitor', np.float32)])
 PC_TYPE = np.dtype([('current', np.float32), ('time', np.int64)])
+LOG_TYPE = np.dtype([('value', np.float32), ('time', np.int64)])
 
 # define the bitmask for individual event filters
 EVENT_BITMASKS = {
@@ -60,6 +61,7 @@ class AmorEventStream:
     packets: np.recarray # PACKET_TYPE
     pulses: np.recarray  # PULSE_TYPE
     proton_current: np.recarray # PC_TYPE
+    device_logs: Dict[str, np.recarray] = field(default_factory=dict) # LOG_TYPE
 
 class EventDatasetProtocol(Protocol):
     """

eos/event_handling.py

@@ -48,7 +48,7 @@ class ApplyParameterOverwrites(EventDataAction):
                     with open(self.config.sampleModel, 'r') as model_yml:
                         model = yaml.safe_load(model_yml)
                 else:
-                    logging.warning(f'  ! the file {self.config.sampleModel}.yml does not exist. Ignored!')
+                    logging.warning(f'  ! the file {self.config.sampleModel} does not exist. Ignored!')
                     return
             else:
                 model = dict(stack=self.config.sampleModel)
@@ -71,55 +71,36 @@ class CorrectSeriesTime(EventDataAction):
         dataset.data.pulses.time -= self.seriesStartTime
         dataset.data.events.wallTime -= self.seriesStartTime
         dataset.data.proton_current.time -= self.seriesStartTime
+        for value in dataset.data.device_logs.values():
+            value.time -= self.seriesStartTime
         start, stop = dataset.data.events.wallTime[0], dataset.data.events.wallTime[-1]
         logging.debug(f'      wall time from {start/1e9:6.1f} s to {stop/1e9:6.1f} s, '
                       f'series time = {self.seriesStartTime/1e9:6.1f}')
 
 
 class AssociatePulseWithMonitor(EventDataAction):
-    def __init__(self, monitorType:MonitorType, lowCurrentThreshold:float):
+    def __init__(self, monitorType:MonitorType):
         self.monitorType = monitorType
-        self.lowCurrentThreshold = lowCurrentThreshold
 
     def perform_action(self, dataset: EventDatasetProtocol)->None:
         logging.debug(f'      using monitor type {self.monitorType}')
         if self.monitorType in [MonitorType.proton_charge or MonitorType.debug]:
-            if not 'wallTime' in dataset.data.events.dtype.names:
-                raise ValueError(
-                    "AssociatePulseWithMonitor requires walltTime to be extracted, please run ExtractWalltime first")
-            monitorPerPulse = self.get_current_per_pulse(dataset.data.pulses.time,
+            dataset.data.pulses.monitor = self.get_current_per_pulse(dataset.data.pulses.time,
                                                               dataset.data.proton_current.time,
                                                               dataset.data.proton_current.current)\
                                                               * 2*dataset.timing.tau * 1e-3
-            # filter low-current pulses
-            dataset.data.pulses.monitor = np.where(
-                    monitorPerPulse > 2*dataset.timing.tau * self.lowCurrentThreshold * 1e-3,
-                    monitorPerPulse, 0)
         elif self.monitorType==MonitorType.time:
             dataset.data.pulses.monitor  = 2*dataset.timing.tau
         else:  # pulses
             dataset.data.pulses.monitor  = 1
 
         if self.monitorType == MonitorType.debug:
+            if not 'wallTime' in dataset.data.events.dtype.names:
+                raise ValueError(
+                    "AssociatePulseWithMonitor requires walltTime for debugging, please run ExtractWalltime first")
             cpp, t_bins = np.histogram(dataset.data.events.wallTime, dataset.data.pulses.time)
             np.savetxt('tme.hst', np.vstack((dataset.data.pulses.time[:-1], cpp, dataset.data.pulses.monitor[:-1])).T)
 
-        if self.monitorType in [MonitorType.proton_charge or MonitorType.debug]:
-            self.monitor_threshold(dataset)
-
-    def monitor_threshold(self, dataset):
-        goodTimeS = dataset.data.pulses.time[dataset.data.pulses.monitor!=0]
-        filter_e = np.logical_not(np.isin(dataset.data.events.wallTime, goodTimeS))
-        dataset.data.events.mask += EVENT_BITMASKS['MonitorThreshold']*filter_e
-        logging.info(f'      low-beam (<{self.lowCurrentThreshold} mC) rejected pulses: '
-                     f'{dataset.data.pulses.monitor.shape[0]-goodTimeS.shape[0]} '
-                     f'out of {dataset.data.pulses.monitor.shape[0]}')
-        logging.info(f'          with {filter_e.sum()} events')
-        if goodTimeS.shape[0]:
-            logging.info(f'      average counts per pulse =  {dataset.data.events.shape[0]/goodTimeS.shape[0]:7.1f}')
-        else:
-            logging.info(f'      average counts per pulse = undefined')
-
     @staticmethod
     def get_current_per_pulse(pulseTimeS, currentTimeS, currents):
         # add currents for early pulses and current time value after last pulse (j+1)
@@ -134,6 +115,28 @@ class AssociatePulseWithMonitor(EventDataAction):
             pulseCurrentS[i] = currents[j]
         return pulseCurrentS
 
+class FilterMonitorThreshold(EventDataAction):
+    def __init__(self, lowCurrentThreshold:float):
+        self.lowCurrentThreshold = lowCurrentThreshold
+
+    def perform_action(self, dataset: EventDatasetProtocol) ->None:
+        if not 'wallTime' in dataset.data.events.dtype.names:
+            raise ValueError(
+                    "FilterMonitorThreshold requires walltTime to be extracted, please run ExtractWalltime first")
+        low_current_filter = dataset.data.pulses.monitor>2*dataset.timing.tau*self.lowCurrentThreshold*1e-3
+        dataset.data.pulses.monitor[np.logical_not(low_current_filter)] = 0.
+        goodTimeS = dataset.data.pulses.time[low_current_filter]
+        filter_e = np.logical_not(np.isin(dataset.data.events.wallTime, goodTimeS))
+
+        dataset.data.events.mask += EVENT_BITMASKS['MonitorThreshold']*filter_e
+        logging.info(f'      low-beam (<{self.lowCurrentThreshold} mC) rejected pulses: '
+                     f'{dataset.data.pulses.monitor.shape[0]-goodTimeS.shape[0]} '
+                     f'out of {dataset.data.pulses.monitor.shape[0]}')
+        logging.info(f'          with {filter_e.sum()} events')
+        if goodTimeS.shape[0]:
+            logging.info(f'      average counts per pulse =  {dataset.data.events.shape[0]/goodTimeS.shape[0]:7.1f}')
+        else:
+            logging.info(f'      average counts per pulse = undefined')
 
 class FilterStrangeTimes(EventDataAction):
     def perform_action(self, dataset: EventDatasetProtocol)->None:
@@ -148,6 +151,9 @@ class TofTimeCorrection(EventDataAction):
         self.correct_chopper_opening = correct_chopper_opening
 
     def perform_action(self, dataset: EventDatasetProtocol) ->None:
+        if not 'delta' in dataset.data.events.dtype.names:
+            raise ValueError(
+                    "TofTimeCorrection requires delta to be extracted, please run AnalyzePixelIDs first")
         d = dataset.data
         if self.correct_chopper_opening:
             d.events.tof -= ( d.events.delta / 180. ) * dataset.timing.tau
@@ -162,7 +168,7 @@ class ApplyMask(EventDataAction):
         # TODO: why is this action time consuming?
         d = dataset.data
         pre_filter = d.events.shape[0]
-        if logging.getLogger().level == logging.DEBUG:
+        if logging.getLogger().level <= logging.DEBUG:
             # only run this calculation if debug level is actually active
             filtered_by_mask = {}
             for key, value in EVENT_BITMASKS.items():
@@ -177,3 +183,20 @@ class ApplyMask(EventDataAction):
             d.events = d.events[fltr]
         post_filter = d.events.shape[0]
         logging.info(f'      number of events: total = {pre_filter:7d}, filtered = {post_filter:7d}')
+        if d.device_logs == {} or not hasattr(dataset, 'update_info_from_logs'):
+            return
+        # filter pulses and logs to allow update of header information
+        from .helpers import add_log_to_pulses
+        times = np.unique(d.events.wallTime)
+        # make sure all log variables are associated with pulses
+        for key, log in d.device_logs.items():
+            if not key in d.pulses.dtype.names:
+                # interpolate the parameter values for all existing pulses
+                add_log_to_pulses(key, dataset)
+        # remove all pulses that have no more events
+        d.pulses = d.pulses[np.isin(d.pulses.time, times)]
+        for key, log in d.device_logs.items():
+            d.device_logs[key] = np.recarray(d.pulses.shape, dtype = log.dtype)
+            d.device_logs[key].time = d.pulses.time
+            d.device_logs[key].value = d.pulses[key]
+        dataset.update_info_from_logs()

eos/file_reader.py

@@ -3,9 +3,9 @@ Reading of Amor NeXus data files to extract metadata and event stream.
 """
 from typing import BinaryIO, List, Union
 
+import sys
 import h5py
 import numpy as np
-import platform
 import logging
 import subprocess
 
@@ -16,7 +16,8 @@ from orsopy.fileio.model_language import SampleModel
 
 from . import const
 from .header import Header
-from .event_data_types import AmorGeometry, AmorTiming, AmorEventStream, PACKET_TYPE, EVENT_TYPE, PULSE_TYPE, PC_TYPE
+from .event_data_types import AmorGeometry, AmorTiming, AmorEventStream, LOG_TYPE, PACKET_TYPE, EVENT_TYPE, PULSE_TYPE, \
+    PC_TYPE
 
 try:
     import zoneinfo
@@ -27,99 +28,145 @@ except ImportError:
 
 # Time zone used to interpret time strings
 AMOR_LOCAL_TIMEZONE = zoneinfo.ZoneInfo(key='Europe/Zurich')
+UTC = zoneinfo.ZoneInfo(key='UTC')
 
-if  platform.node().startswith('amor'):
-    NICOS_CACHE_DIR = '/home/amor/nicosdata/amor/cache/'
-    GREP = '/usr/bin/grep "%s"'
-else:
-    NICOS_CACHE_DIR = None
-
-class AmorEventData:
+class AmorHeader:
     """
-    Read one amor NeXus datafile and extract relevant header information.
-
-    Implements EventDatasetProtocol
+    Collects header information from Amor NeXus fiel without reading event data.
     """
-    file_list: List[str]
-    first_index: int
-    last_index: int = -1
-    EOF: bool = False
-    max_events: int
-    owner: fileio.Person
-    experiment: fileio.Experiment
-    sample: fileio.Sample
-    instrument_settings: fileio.InstrumentSettings
-    geometry: AmorGeometry
-    timing: AmorTiming
-    data: AmorEventStream
+    # mapping of names to (hdf_path, dtype, nicos_key[, suffix])
+    hdf_paths = dict(
+            title=('entry1/title', str),
+            proposal_id=('entry1/proposal_id', str),
+            user_name=('entry1/user/name', str),
+            user_email=('entry1/user/email', str),
+            sample_name=('entry1/sample/name', str),
+            source_name=('entry1/Amor/source/name', str),
+            sample_model=('entry1/sample/model', str),
+            start_time=('entry1/start_time', str),
+            start_time_fallback=('entry1/Amor/instrument_control_parameters/start_time', str),
 
-    eventStartTime: np.int64
+            chopper_separation=('entry1/Amor/chopper/pair_separation', float),
+            detector_distance=('entry1/Amor/detector/transformation/distance', float),
+            chopper_distance=('entry1/Amor/chopper/distance', float),
+            sample_temperature=('entry1/sample/temperature', float),
+            sample_magnetic_field=('entry1/sample/magnetic_field', float),
 
-    def __init__(self, fileName:Union[str, h5py.File, BinaryIO], first_index:int=0, max_events:int=100_000_000):
+            mu=('entry1/Amor/instrument_control_parameters/mu', float, 'mu'),
+            nu=('entry1/Amor/instrument_control_parameters/nu', float, 'nu'),
+            kap=('entry1/Amor/instrument_control_parameters/kappa', float, 'kappa'),
+            kad=('entry1/Amor/instrument_control_parameters/kappa_offset', float, 'kappa_offset'),
+            div=('entry1/Amor/instrument_control_parameters/div', float, 'div'),
+            ch1_trigger_phase=('entry1/Amor/chopper/ch1_trigger_phase', float, 'ch1_trigger_phase'),
+            ch2_trigger_phase=('entry1/Amor/chopper/ch2_trigger_phase', float, 'ch2_trigger_phase'),
+            chopper_speed=('entry1/Amor/chopper/rotation_speed', float, 'chopper_phase'),
+            chopper_phase=('entry1/Amor/chopper/phase', float, 'chopper_phase'),
+            polarization_config_label=('entry1/Amor/polarization/configuration', int, 'polarization_config_label', '/*'),
+            )
+
+    def __init__(self, fileName:Union[str, h5py.File, BinaryIO]):
         if type(fileName) is str:
-            self.file_list = [fileName]
+            logging.warning(f'    {fileName.split("/")[-1]}')
             self.hdf = h5py.File(fileName, 'r', swmr=True)
         elif type(fileName) is h5py.File:
-            self.file_list = [fileName.filename]
             self.hdf = fileName
         else:
-            self.file_list = [repr(fileName)]
             self.hdf = h5py.File(fileName, 'r')
-        self.first_index = first_index
-        self.max_events = max_events
+
+        self._log_keys = []
 
         self.read_header_info()
         self.read_instrument_configuration()
-        self.read_event_stream()
-
-        # actions applied to any dataset
-        self.read_chopper_trigger_stream()
-        self.read_proton_current_stream()
 
         if type(fileName) is str:
             # close the input file to free memory, only if the file was opened in this object
             self.hdf.close()
         del(self.hdf)
 
-    def _replace_if_missing(self, key, nicos_key, dtype=float):
+    def _replace_if_missing(self, key, nicos_key, dtype=float, suffix=''):
+        from .nicos_interface import lookup_nicos_value
+        year = self.fileDate.strftime('%Y')
+        return lookup_nicos_value(key, nicos_key, dtype, suffix, year)
+
+    def rv(self, key):
+        """
+        Generic read value methos based on key in hdf_paths dictionary.
+        """
+        hdf_path, dtype, *nicos = self.hdf_paths[key]
         try:
-            return dtype(self.hdf[f'/entry1/Amor/{key}'][0])
-        except(KeyError, IndexError):
-            if NICOS_CACHE_DIR:
+            hdfgrp = self.hdf[hdf_path]
+            if hdfgrp.attrs.get('NX_class', None) == 'NXlog':
+                # use the last value that was recoreded before the count started
+                time_column = hdfgrp['time'][:]
                 try:
-                    logging.warning(f"     using parameter {nicos_key} from nicos cache")
-                    year_date = self.fileDate.strftime('%Y')
-                    value = str(subprocess.getoutput(f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year_date}')).split('\t')[-1]
-                    return dtype(value)
-                except Exception:
-                    logging.error("Couldn't get value from nicos cache", exc_info=True)
-                    return dtype(0)
+                    start_index = np.where(time_column<=self._start_time_ns)[0][0]
+                except IndexError:
+                    start_index = 0
+                if hdfgrp['value'].ndim==1:
+                    output = dtype(hdfgrp['value'][start_index])
+                else:
+                    output = dtype(hdfgrp['value'][start_index, 0])
+                # make sure key is only appended if no exception was raised
+                self._log_keys.append(key)
+                return output
+            elif dtype is str:
+                return self.read_string(hdf_path)
             else:
-                logging.warning(f"     parameter {key} not found, relpace by zero")
-                return dtype(0)
+                if len(hdfgrp.shape)==1:
+                    return dtype(hdfgrp[0])
+                else:
+                    return dtype(hdfgrp[()])
+        except (KeyError, IndexError):
+            if nicos:
+                nicos_key = nicos[0]
+                suffix = nicos[1] if len(nicos)>1 else ''
+                return self._replace_if_missing(key, nicos_key, dtype, suffix)
+            else:
+                raise
+
+
+    def read_string(self, path):
+        if not np.shape(self.hdf[path]):
+            return self.hdf[path][()].decode('utf-8')
+        else:
+            # format until 2025
+            return self.hdf[path][0].decode('utf-8')
 
     def read_header_info(self):
+        self._start_time_ns = np.uint64(0)
+        try:
+            start_time = self.rv('start_time')
+        except KeyError:
+            start_time = self.rv('start_time_fallback')
+
+        # extract start time as unix time, adding UTC offset of 1h to time string
+        if start_time.endswith('Z') and sys.version_info.minor<11:
+            # older python versions did not support Z format
+            start_time = start_time[:-1]
+            TZ = UTC
+        else:
+            TZ = AMOR_LOCAL_TIMEZONE
+        start_date = datetime.fromisoformat(start_time)
+        self.fileDate = start_date.replace(tzinfo=TZ)
+        self._start_time_ns = np.uint64(self.fileDate.timestamp()*1e9)
+
         # 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.rv('title')
+        proposal_id = self.rv('proposal_id')
+        user_name = self.rv('user_name')
         user_affiliation = 'unknown'
-        user_email = self.hdf['entry1/user/email'][0].decode('utf-8')
+        user_email = self.rv('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.rv('sample_name')
+        instrumentName = 'Amor'
+        source = self.rv('source_name')
+        sourceProbe = 'neutron'
+        model = self.rv('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')
-        # 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)
 
         self.owner = fileio.Person(
                 name=user_name,
@@ -137,26 +184,42 @@ class AmorEventData:
                 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:
+            try:
+                sample_temperature = self.rv('sample_temperature')
+            except IndexError: pass
+            else:
+                self.sample.sample_parameters['temperature'] = fileio.Value(sample_temperature, unit='K')
+        if self.hdf['entry1/sample'].get('magnetic_field', None) is not None:
+            try:
+                sample_magnetic_field = self.rv('sample_magnetic_field')
+            except IndexError: pass
+            else:
+                self.sample.sample_parameters['magnetic_field'] = fileio.Value(sample_magnetic_field, unit='T')
 
     def read_instrument_configuration(self):
-        chopperSeparation = float(np.take(self.hdf['entry1/Amor/chopper/pair_separation'], 0))
-        detectorDistance = float(np.take(self.hdf['entry1/Amor/detector/transformation/distance'], 0))
-        chopperDetectorDistance = detectorDistance-float(np.take(self.hdf['entry1/Amor/chopper/distance'], 0))
+        chopperSeparation = self.rv('chopper_separation')
+        detectorDistance = self.rv('detector_distance')
+        chopperDistance = self.rv('chopper_distance')
+        chopperDetectorDistance = detectorDistance - chopperDistance
 
-        polarizationConfigs = ['unpolarized', 'unpolarized', 'po', 'mo', 'op', 'pp', 'mp', 'om', 'pm', 'mm']
-
-        mu = self._replace_if_missing('instrument_control_parameters/mu', 'mu', float)
-        nu = self._replace_if_missing('instrument_control_parameters/nu', 'nu', float)
-        kap = self._replace_if_missing('instrument_control_parameters/kappa', 'kappa', float)
-        kad = self._replace_if_missing('instrument_control_parameters/kappa_offset', 'kad', float)
-        div = self._replace_if_missing('instrument_control_parameters/div', 'div', float)
-        ch1TriggerPhase = self._replace_if_missing('chopper/ch1_trigger_phase', 'ch1_trigger_phase', float)
-        ch2TriggerPhase = self._replace_if_missing('chopper/ch2_trigger_phase', 'ch2_trigger_phase', float)
+        mu = self.rv('mu')
+        nu = self.rv('nu')
+        kap = self.rv('kap')
+        kad = self.rv('kad')
+        div = self.rv('div')
+        ch1TriggerPhase = self.rv('ch1_trigger_phase')
+        ch2TriggerPhase = self.rv('ch2_trigger_phase')
         try:
             chopperTriggerTime = (float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][7]) \
                                   -float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][0])) \
@@ -164,8 +227,8 @@ class AmorEventData:
             chopperTriggerTimeDiff = float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][2])
         except (KeyError, IndexError):
             logging.debug('      chopper speed and phase taken from .hdf file')
-            chopperSpeed = self._replace_if_missing('chopper/rotation_speed', 'chopper_phase', float)
-            chopperPhase = self._replace_if_missing('chopper/phase', 'chopper_phase', float)
+            chopperSpeed = self.rv('chopper_speed')
+            chopperPhase = self.rv('chopper_phase')
             tau = 30/chopperSpeed
         else:
             tau = int(1e-6*chopperTriggerTime/2+0.5)*(1e-3)
@@ -177,8 +240,8 @@ class AmorEventData:
                                      chopperSeparation, detectorDistance, chopperDetectorDistance)
         self.timing = AmorTiming(ch1TriggerPhase, ch2TriggerPhase, chopperSpeed, chopperPhase, tau)
 
-        polarizationConfigLabel = self._replace_if_missing('polarization/configuration/average_value', 'polarizer_config_label', int)
-        polarizationConfig = fileio.Polarization(polarizationConfigs[polarizationConfigLabel])
+        polarizationConfigLabel = self.rv('polarization_config_label')
+        polarizationConfig = fileio.Polarization(const.polarizationConfigs[polarizationConfigLabel])
         logging.debug(f'      polarization configuration: {polarizationConfig} (index {polarizationConfigLabel})')
 
 
@@ -187,9 +250,11 @@ class AmorEventData:
                                                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,6 +288,58 @@ class AmorEventData:
         header.sample = self.sample
         header.measurement_instrument_settings = self.instrument_settings
 
+
+class AmorEventData(AmorHeader):
+    """
+    Read one amor NeXus datafile and extract relevant header information.
+
+    Implements EventDatasetProtocol
+    """
+    file_list: List[str]
+    first_index: int
+    last_index: int = -1
+    EOF: bool = False
+    max_events: int
+    owner: fileio.Person
+    experiment: fileio.Experiment
+    sample: fileio.Sample
+    instrument_settings: fileio.InstrumentSettings
+    geometry: AmorGeometry
+    timing: AmorTiming
+    data: AmorEventStream
+
+    eventStartTime: np.int64
+
+    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'    {fileName.split("/")[-1]}')
+            self.file_list = [fileName]
+            hdf = h5py.File(fileName, 'r', swmr=True)
+        elif type(fileName) is h5py.File:
+            self.file_list = [fileName.filename]
+            hdf = fileName
+        else:
+            self.file_list = [repr(fileName)]
+            hdf = h5py.File(fileName, 'r')
+        self.first_index = first_index
+        self.max_events = max_events
+
+        super().__init__(hdf)
+        self.hdf = hdf
+        try:
+            self.read_event_stream()
+        except EOFError:
+            self.hdf.close()
+            del(self.hdf)
+            raise
+        self.read_log_streams()
+
+        if type(fileName) is str:
+            # close the input file to free memory, only if the file was opened in this object
+            self.hdf.close()
+        del(self.hdf)
+
+
     def read_event_stream(self):
         """
         Read the actual event data from file. If file is too large, find event index from packets
@@ -230,22 +347,34 @@ class AmorEventData:
         """
         packets = np.recarray(self.hdf['/entry1/Amor/detector/data/event_index'].shape, dtype=PACKET_TYPE)
         packets.start_index = self.hdf['/entry1/Amor/detector/data/event_index'][:]
-        packets.Time = self.hdf['/entry1/Amor/detector/data/event_time_zero'][:]
+        packets.time = self.hdf['/entry1/Amor/detector/data/event_time_zero'][:]
         try:
             # packet index that matches first event index
             start_packet = int(np.where(packets.start_index==self.first_index)[0][0])
         except IndexError:
-            raise IndexError(f'No event packet found starting at event #{self.first_index}')
+            raise EOFError(f'No event packet found starting at event #{self.first_index}, '
+                           f'number of events is {self.hdf["/entry1/Amor/detector/data/event_time_offset"].shape[0]}')
         packets = packets[start_packet:]
+        if packets.shape[0]==0:
+            raise EOFError(f'No more packets left after start_packet filter')
 
         nevts = self.hdf['/entry1/Amor/detector/data/event_time_offset'].shape[0]
         if (nevts-self.first_index)>self.max_events:
             end_packet = np.where(packets.start_index<=(self.first_index+self.max_events))[0][-1]
-            self.last_index = packets.start_index[end_packet]-1
+            end_packet = max(1, end_packet)
+            if len(packets)==1:
+                self.last_index = nevts-1
+            else:
+                self.last_index = packets.start_index[end_packet]-1
             packets = packets[:end_packet]
         else:
             self.last_index = nevts-1
             self.EOF = True
+
+        if packets.shape[0]==0:
+            raise EOFError(f'No more packets left after end_packet filter, first_index={self.first_index}, '
+                           f'max_events={self.max_events}, nevts={nevts}')
+
         nevts = self.last_index+1-self.first_index
 
         # adapte packet to event index relation
@@ -256,15 +385,54 @@ class AmorEventData:
         events.pixelID = self.hdf['/entry1/Amor/detector/data/event_id'][self.first_index:self.last_index+1]
         events.mask = 0
 
-        pulses = self.read_chopper_trigger_stream()
-        current = self.read_proton_current_stream()
+        pulses = self.read_chopper_trigger_stream(packets)
+        current = self.read_proton_current_stream(packets)
         self.data = AmorEventStream(events, packets, pulses, current)
 
         if self.first_index>0 and not self.EOF:
             # label the file name if not all events were used
             self.file_list[0] += f'[{self.first_index}:{self.last_index}]'
 
-    def read_chopper_trigger_stream(self):
+    def read_log_streams(self):
+        """
+        Read the individual NXlog datasets that can later be used for filtering etc.
+        """
+        for key in self._log_keys:
+            hdf_path, dtype, *_ = self.hdf_paths[key]
+            hdfgroup = self.hdf[hdf_path]
+            shape = hdfgroup['time'].shape
+            data = np.recarray(shape, dtype=np.dtype([('value', self.hdf_paths[key][1]), ('time', np.int64)]))
+            data.time = hdfgroup['time'][:]
+            if len(hdfgroup['value'].shape)==1:
+                data.value = hdfgroup['value'][:]
+            else:
+                data.value = hdfgroup['value'][:, 0]
+            self.data.device_logs[key] = data
+
+    def update_info_from_logs(self):
+        RELEVANT_ITEMS = ['sample_temperature', 'sample_magnetic_field', 'polarization_config_label']
+        for key, log in self.data.device_logs.items():
+            if key not in RELEVANT_ITEMS:
+                continue
+            if log.value.dtype in [np.int8, np.int16, np.int32, np.int64]:
+                # for integer items (flags) report the most common one
+                value = np.bincount(log.value).argmax()
+                if logging.getLogger().getEffectiveLevel() <= logging.DEBUG \
+                        and np.unique(log.value).shape[0]>1:
+                    logging.debug(f'        filtered values for {key} not unique, '
+                                  f'has {np.unique(log.value).shape[0]} values')
+            else:
+                value = log.value.mean()
+            if key == 'polarization_config_label':
+                self.instrument_settings.polarization = fileio.Polarization(const.polarizationConfigs[value])
+            elif key == 'sample_temperature':
+                self.sample.sample_parameters['temperature'].magnitue = value
+            elif key == 'sample_magnetic_field':
+                self.sample.sample_parameters['magnetic_field'].magnitue = value
+
+
+
+    def read_chopper_trigger_stream(self, packets):
         chopper1TriggerTime = np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][:-2], dtype=np.int64)
         #self.chopper2TriggerTime = self.chopper1TriggerTime + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time'][:-2], dtype=np.int64)
         #                           + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][:], dtype=np.int64)
@@ -272,7 +440,7 @@ class AmorEventData:
             startTime = chopper1TriggerTime[0]
             pulseTimeS = chopper1TriggerTime
         else:
-            logging.warn('     no chopper trigger data available, using event steram instead')
+            logging.critical('     No chopper trigger data available, using event steram instead, pulse filtering will fail!')
             startTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][0], dtype=np.int64)
             stopTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-2], dtype=np.int64)
             pulseTimeS = np.arange(startTime, stopTime, self.timing.tau*1e9, dtype=np.int64)
@@ -280,18 +448,22 @@ class AmorEventData:
         pulses.time = pulseTimeS
         pulses.monitor = 1. # default is monitor pulses as it requires no calculation
         # apply filter in case the events were filtered
-        if self.first_index>0 or not self.EOF:
-            pulses = pulses[(pulses.time>=self.data.packets.Time[0])&(pulses.time<=self.data.packets.Time[-1])]
+        if (self.first_index>0 or not self.EOF):
+            pulses = pulses[(pulses.time>=packets.time[0])&(pulses.time<=packets.time[-1])]
         self.eventStartTime = startTime
         return pulses
 
-    def read_proton_current_stream(self):
+    def read_proton_current_stream(self, packets):
         proton_current = np.recarray(self.hdf['entry1/Amor/detector/proton_current/time'].shape, dtype=PC_TYPE)
         proton_current.time = self.hdf['entry1/Amor/detector/proton_current/time'][:]
-        proton_current.current = self.hdf['entry1/Amor/detector/proton_current/value'][:,0]
+        if self.hdf['entry1/Amor/detector/proton_current/value'].ndim==1:
+            proton_current.current = self.hdf['entry1/Amor/detector/proton_current/value'][:]
+        else:
+            proton_current.current = self.hdf['entry1/Amor/detector/proton_current/value'][:,0]
+
         if self.first_index>0 or not self.EOF:
-            proton_current = proton_current[(proton_current.time>=self.data.packets.Time[0])&
-                                            (proton_current.time<=self.data.packets.Time[-1])]
+            proton_current = proton_current[(proton_current.time>=packets.time[0])&
+                                            (proton_current.time<=packets.time[-1])]
         return proton_current
 
     def info(self):

eos/header.py

@@ -5,6 +5,7 @@ Class to handle Orso header information that changes gradually during the reduct
 import platform
 import sys
 from datetime import datetime
+from typing import List, Literal
 
 from orsopy import fileio
 
@@ -13,13 +14,16 @@ from . import __version__
 
 class Header:
     """orso compatible output file header content"""
+    owner: fileio.Person
+    experiment: fileio.Experiment
+    sample: fileio.Sample
+    measurement_instrument_settings: fileio.InstrumentSettings
+    measurement_scheme: Literal["angle- and energy-dispersive", "angle-dispersive", "energy-dispersive"]
+    measurement_data_files: List[fileio.File]
+    measurement_additional_files: List[fileio.File]
+
 
     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    = []
 
@@ -64,10 +68,3 @@ class Header:
         if columns is None:
             columns = self.columns()
         return fileio.Orso(ds, red, columns+extra_columns)
-    #-------------------------------------------------------------------------------------------------
-    def create_call_string(self):
-        callString = ' '.join(sys.argv)
-        if '-Y' not in callString:
-            callString += f' -Y {datetime.now().year}'
-        return callString
-    #-------------------------------------------------------------------------------------------------

eos/helpers.py

@@ -1,10 +1,35 @@
-"""
-Helper functions used during calculations. Uses numba enhanced functions if available, otherwise numpy based
-fallback is imported.
-"""
-
-try:
-    from .helpers_numba import merge_frames, extract_walltime, filter_project_x, calculate_derived_properties_focussing
-except ImportError:
-    from .helpers_fallback import merge_frames, extract_walltime, filter_project_x, calculate_derived_properties_focussing
-
+"""
+Helper functions used during calculations. Uses numba enhanced functions if available, otherwise numpy based
+fallback is imported.
+"""
+import numpy as np
+from .event_data_types import EventDatasetProtocol, append_fields
+
+try:
+    from .helpers_numba import merge_frames, extract_walltime, filter_project_x, calculate_derived_properties_focussing
+except ImportError:
+    from .helpers_fallback import merge_frames, extract_walltime, filter_project_x, calculate_derived_properties_focussing
+
+def add_log_to_pulses(key, dataset: EventDatasetProtocol):
+    """
+    Add a log value for each pulse to the pulses array.
+    """
+    pulses = dataset.data.pulses
+    log_data = dataset.data.device_logs[key]
+    if log_data.time[0]>0:
+        logTimeS = np.hstack([[0], log_data.time, [pulses.time[-1]+1]])
+        logValues = np.hstack([[log_data.value[0]], log_data.value])
+    else:
+        logTimeS = np.hstack([log_data.time, [pulses.time[-1]+1]])
+        logValues = log_data.value
+    pulseLogS = np.zeros(pulses.time.shape[0], dtype=log_data.value.dtype)
+    j = 0
+    for i, ti in enumerate(pulses.time):
+        # find the last current item that was before this pulse
+        while ti>=logTimeS[j+1]:
+            j += 1
+        pulseLogS[i] = logValues[j]
+    pulses = append_fields(pulses, [(key, pulseLogS.dtype)])
+    pulses[key] = pulseLogS
+    dataset.data.pulses = pulses
+

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,13 +97,14 @@ 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
 
     @cache
     def z(self):
+        # TODO: shouldn't this be -0.5 to be the edges of each pixel?
         return np.arange(Detector.nBlades*Detector.nWires+1)
 
     @cache

eos/kafka_events.py

@@ -0,0 +1,165 @@
+"""
+Collect AMOR detector events send via Kafka.
+"""
+
+import logging
+import numpy as np
+from threading import Thread, Event
+from time import time
+
+from .event_data_types import AmorGeometry, AmorTiming, AmorEventStream, PACKET_TYPE, EVENT_TYPE, PULSE_TYPE, PC_TYPE
+
+from uuid import uuid4
+from streaming_data_types.eventdata_ev44 import EventData
+from streaming_data_types.logdata_f144 import ExtractedLogData
+from streaming_data_types import deserialise_f144, deserialise_ev44
+from confluent_kafka import Consumer
+
+from .header import Header
+
+
+try:
+    from streaming_data_types.utils import get_schema
+except ImportError:
+    from streaming_data_types.utils import _get_schema as get_schema
+
+
+KAFKA_BROKER = 'linkafka01.psi.ch:9092'
+AMOR_EVENTS = 'amor_detector'
+AMOR_NICOS = 'amor_nicosForwarder'
+
+class KafkaFrozenData:
+    """
+    Represents event stream data from Kafka at a given time.
+    Will be returned by KafkaEventData to be use in conjunction
+    with data processing and projections.
+
+    Implements EventDatasetProtocol
+    """
+    geometry: AmorGeometry
+    timing: AmorTiming
+    data: AmorEventStream
+
+    def __init__(self, geometry, timing, data, monitor=1.):
+        self.geometry = geometry
+        self.timing = timing
+        self.data = data
+        self.monitor = monitor
+
+    def append(self, other):
+        raise NotImplementedError("can't append live datastream to other event data")
+
+    def update_header(self, header:Header):
+        # maybe makes sense later, but for now just used for live vizualization
+        ...
+
+class KafkaEventData(Thread):
+    """
+    Read Nicos information and events from Kafka. Creates a background
+    thread that listens to Kafka events and converts them to eos compatible information.
+    """
+    geometry: AmorGeometry
+    timing: AmorTiming
+    events: np.recarray
+
+    def __init__(self):
+        self.stop_event = Event()
+        self.stop_counting = Event()
+        self.new_events = Event()
+        self.last_read = 0
+        self.last_read_time = 0.
+        self.start_time = time()
+        self.consumer = Consumer(
+                {'bootstrap.servers': 'linkafka01.psi.ch:9092',
+                 'group.id': uuid4()})
+        self.consumer.subscribe([AMOR_EVENTS, AMOR_NICOS])
+        self.geometry = AmorGeometry(1.0, 2.0, 0., 0., 1.5, 10.0, 4.0, 10.0)
+        self.timing = AmorTiming(0., 0., 500., 0., 30./500.)
+        # create empty dataset
+        self.events = np.recarray(0, dtype=EVENT_TYPE)
+        super().__init__()
+
+    def run(self):
+        while not self.stop_event.is_set():
+            messages = self.consumer.consume(10, timeout=1)
+            for message in messages:
+                self.process_message(message)
+
+    def process_message(self, message):
+        if message.error():
+            logging.info(f"  received Kafka message with error: {message.error()}")
+            return
+        schema = get_schema(message.value())
+        if message.topic()==AMOR_EVENTS and schema=='ev44':
+            events:EventData = deserialise_ev44(message.value())
+            self.add_events(events)
+            self.new_events.set()
+            logging.debug(f'  new events {events}')
+        elif message.topic()==AMOR_NICOS and schema=='f144':
+            nicos_data:ExtractedLogData = deserialise_f144(message.value())
+            if nicos_data.source_name in self.nicos_mapping.keys():
+                logging.debug(f'  {nicos_data.source_name} = {nicos_data.value}')
+                self.update_instrument(nicos_data)
+
+    def add_events(self, events:EventData):
+        """
+        Add new events to the Dataset. The object keeps raw events
+        and only copies the latest set to the self.data object,
+        this allows to run the event processing to be performed on a "clean"
+        evnet stream each time.
+        """
+        if self.stop_counting.is_set():
+            return
+        prev_size = self.events.shape[0]
+        new_events = events.pixel_id.shape[0]
+        self.events.resize(prev_size+new_events, refcheck=False)
+        self.events.pixelID[prev_size:] = events.pixel_id
+        self.events.mask[prev_size:] = 0
+        self.events.tof[prev_size:] = events.time_of_flight/1.e9
+
+    nicos_mapping = {
+        'mu': ('geometry', 'mu'),
+        'nu':  ('geometry', 'nu'),
+        'kappa':  ('geometry', 'kap'),
+        'kappa_offset':  ('geometry', 'kad'),
+        'ch1_trigger_phase':  ('timing', 'ch1TriggerPhase'),
+        'ch2_trigger_phase':  ('timing', 'ch2TriggerPhase'),
+        'ch2_speed':  ('timing', 'chopperSpeed'),
+        'chopper_phase':  ('timing', 'chopperPhase'),
+        }
+
+    def update_instrument(self, nicos_data:ExtractedLogData):
+        if nicos_data.source_name in self.nicos_mapping:
+            attr, subattr = self.nicos_mapping[nicos_data.source_name]
+            setattr(getattr(self, attr), subattr, nicos_data.value)
+            if nicos_data.source_name=='ch2_speed':
+                self.timing.tau = 30./self.timing.chopperSpeed
+
+    def monitor(self):
+        return time()-self.start_time
+
+    def restart(self):
+        # empty event buffer
+        self.events = np.recarray(0, dtype=EVENT_TYPE)
+        self.stop_counting.clear()
+        self.last_read = 0
+        self.start_time = time()
+        self.new_events.clear()
+
+    def get_events(self, total_counts=False):
+        packets = np.recarray(0, dtype=PACKET_TYPE)
+        pulses = np.recarray(0, dtype=PULSE_TYPE)
+        pc = np.recarray(0, dtype=PC_TYPE)
+        if total_counts:
+            last_read = 0
+        else:
+            last_read = self.last_read
+        if last_read>=self.events.shape[0]:
+            raise EOFError("No new events arrived")
+        data = AmorEventStream(self.events[last_read:].copy(), packets, pulses, pc)
+        self.last_read = self.events.shape[0]
+        self.new_events.clear()
+        t_now = time()
+        monitor = t_now-self.last_read_time
+        self.last_read_time = t_now
+        return KafkaFrozenData(self.geometry, self.timing, data, monitor=monitor)

eos/kafka_serializer.py

@@ -0,0 +1,283 @@
+"""
+Allows to send eos projections to Kafka using ESS histogram serialization.
+
+For histogram_h01 the message is build using:
+
+hist = {
+    "source": "some_source",
+    "timestamp": 123456,
+    "current_shape": [2, 5],
+    "dim_metadata": [
+        {
+            "length": 2,
+            "unit": "a",
+            "label": "x",
+            "bin_boundaries": np.array([10, 11, 12]),
+        },
+        {
+            "length": 5,
+            "unit": "b",
+            "label": "y",
+            "bin_boundaries": np.array([0, 1, 2, 3, 4, 5]),
+        },
+    ],
+    "last_metadata_timestamp": 123456,
+    "data": np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]]),
+    "errors": np.array([[5, 4, 3, 2, 1], [10, 9, 8, 7, 6]]),
+    "info": "info_string",
+}
+"""
+import logging
+from typing import List, Tuple, Union
+from threading import Thread, Event
+
+import numpy as np
+import json
+from time import time
+from dataclasses import dataclass, asdict
+from streaming_data_types import histogram_hs01
+from confluent_kafka import Producer, Consumer
+
+from uuid import uuid4
+
+from .projection import TofZProjection, YZProjection
+
+KAFKA_BROKER = 'linkafka01.psi.ch:9092'
+KAFKA_TOPICS = {
+    'histogram': 'amor_histograms',
+    'response': 'amor_histResponse',
+    'command':  'amor_histCommands'
+    }
+
+def ktime():
+    return int(time()*1_000)
+
+@dataclass
+class DimMetadata:
+    length: int
+    unit: str
+    label: str
+    bin_boundaries: np.ndarray
+
+@dataclass
+class HistogramMessage:
+    source: str
+    timestamp: int
+    current_shape: Tuple[int, int]
+    dim_metadata: Tuple[DimMetadata, DimMetadata]
+    last_metadata_timestamp: int
+    data: np.ndarray
+    errors: np.ndarray
+    info: str
+
+    def serialize(self):
+        return histogram_hs01.serialise_hs01(asdict(self))
+
+@dataclass
+class CommandMessage:
+    msg_id: str
+
+    cmd=None
+
+    @classmethod
+    def get_message(cls, data):
+        """
+        Uses the sub-class cmd attribute to select which message to retugn
+        """
+        msg = dict([(ci.cmd, ci) for ci in cls.__subclasses__()])
+        return msg[data['cmd']](**data)
+
+
+@dataclass
+class Stop(CommandMessage):
+    hist_id: str
+    id: str
+    cmd:str = 'stop'
+
+@dataclass
+class HistogramConfig:
+    id: str
+    type: str
+    data_brokers: List[str]
+    topic: str
+    data_topics: List[str]
+    tof_range: Tuple[float, float]
+    det_range: Tuple[int, int]
+    num_bins: int
+    width: int
+    height: int
+    left_edges: list
+    source: str
+
+@dataclass
+class ConfigureHistogram(CommandMessage):
+    histograms: List[HistogramConfig]
+    start: int
+    cmd:str = 'config'
+
+    def __post_init__(self):
+        self.histograms = [HistogramConfig(**cfg) for cfg in self.histograms]
+
+
+class ESSSerializer:
+
+    def __init__(self):
+        self.producer = Producer({
+            'bootstrap.servers': KAFKA_BROKER,
+            'message.max.bytes': 4_000_000,
+            })
+        self.consumer = Consumer({
+            'bootstrap.servers': KAFKA_BROKER,
+            "group.id": uuid4(),
+            "default.topic.config": {"auto.offset.reset": "latest"},
+            })
+        self._active_histogram_yz = None
+        self._active_histogram_tofz = None
+        self.new_count_started = Event()
+        self.count_stopped = Event()
+
+        self.consumer.subscribe([KAFKA_TOPICS['command']])
+
+    def process_message(self, message):
+        if message.error():
+            logging.error("Command Consumer Error: %s", message.error())
+        else:
+            command = json.loads(message.value().decode())
+            try:
+                command = CommandMessage.get_message(command)
+            except Exception:
+                logging.error(f'Could not interpret message: \n{command}', exc_info=True)
+                return
+            logging.info(command)
+            resp = json.dumps({
+                "msg_id":   getattr(command, "id", None) or command.msg_id,
+                "response": "ACK",
+                "message":  ""
+                })
+            self.producer.produce(
+                    topic=KAFKA_TOPICS['response'],
+                    value=resp
+                    )
+            self.producer.flush()
+            if isinstance(command, Stop):
+                if command.hist_id in [self._active_histogram_yz, self._active_histogram_tofz]:
+                    self.count_stopped.set()
+                else:
+                    return
+            elif isinstance(command, ConfigureHistogram):
+                for hist in command.histograms:
+                    if hist.topic == KAFKA_TOPICS['histogram']+'_YZ':
+                        self._active_histogram_yz = hist.id
+                        logging.debug(f"   histogram data_topic: {hist.data_topics}")
+                        self._start = command.start
+                        self.count_stopped.clear()
+                        self.new_count_started.set()
+                    if hist.topic == KAFKA_TOPICS['histogram']+'_TofZ':
+                        self._active_histogram_tofz = hist.id
+
+    def receive(self, timeout=5):
+        rec = self.consumer.poll(timeout)
+        if rec is not None:
+            self.process_message(rec)
+            return True
+        else:
+            return False
+
+    def receive_loop(self):
+        while not self._stop_receiving.is_set():
+            try:
+                self.receive()
+            except Exception:
+                logging.error("Exception while receiving", exc_info=True)
+
+    def start_command_thread(self):
+        self._stop_receiving = Event()
+        self._command_thread = Thread(target=self.receive_loop)
+        self._command_thread.start()
+
+    def end_command_thread(self, event=None):
+        self._stop_receiving.set()
+        self._command_thread.join()
+
+    def acked(self, err, msg):
+        # We need to have callback to produce-method to catch server errors
+        if err is not None:
+            logging.warning("Failed to deliver message: %s: %s" % (str(msg), str(err)))
+        else:
+            logging.debug("Message produced: %s" % (str(msg)))
+
+    def send(self, proj: Union[YZProjection, TofZProjection], final=False):
+        if final:
+            state = 'FINISHED'
+        else:
+            state = 'COUNTING'
+        if isinstance(proj, YZProjection):
+            if self._active_histogram_yz is None:
+                return
+            suffix = 'YZ'
+            message = HistogramMessage(
+                source='amor-eos',
+                timestamp=ktime(),
+                current_shape=(proj.y.shape[0]-1, proj.z.shape[0]-1),
+                dim_metadata=(
+                    DimMetadata(
+                            length=proj.y.shape[0]-1,
+                            unit="pixel",
+                            label="Y",
+                            bin_boundaries=proj.y,
+                            ),
+                    DimMetadata(
+                            length=proj.z.shape[0]-1,
+                            unit="pixel",
+                            label="Z",
+                            bin_boundaries=proj.z,
+                            )
+                ),
+                last_metadata_timestamp=0,
+                data=proj.data.cts,
+                errors=np.sqrt(proj.data.cts),
+                info=json.dumps({
+                    "start": self._start,
+                    "state": state,
+                    "num events": proj.data.cts.sum()
+                })
+                )
+            logging.info(f"   {state}: Sending {proj.data.cts.sum()} events to Nicos")
+        elif isinstance(proj, TofZProjection):
+            if self._active_histogram_tofz is None:
+                return
+            suffix = 'TofZ'
+            message = HistogramMessage(
+                source='amor-eos',
+                timestamp=ktime(),
+                current_shape=(proj.tof.shape[0]-1, proj.z.shape[0]-1),
+                dim_metadata=(
+                    DimMetadata(
+                            length=proj.tof.shape[0]-1,
+                            unit="ms",
+                            label="ToF",
+                            bin_boundaries=proj.tof,
+                            ),
+                    DimMetadata(
+                            length=proj.z.shape[0]-1,
+                            unit="pixel",
+                            label="Z",
+                            bin_boundaries=proj.z,
+                            ),
+                ),
+                last_metadata_timestamp=0,
+                data=proj.data.cts,
+                errors=np.sqrt(proj.data.cts),
+                info=json.dumps({
+                    "start": self._start,
+                    "state": state,
+                    "num events": proj.data.I.sum()
+                })
+                )
+        else:
+            raise NotImplementedError(f"Histogram for {proj.__class__.__name__} not implemented")
+
+        self.producer.produce(value=message.serialize(),
+                              topic=KAFKA_TOPICS['histogram']+'_'+suffix,
+                              callback=self.acked)
+        self.producer.flush()

eos/ls.py

@@ -0,0 +1,54 @@
+"""
+eosls executable script to list available datafiles in current folder with some metadata information.
+
+Author: Jochen Stahn (algorithms, python draft),
+        Artur Glavic (structuring and optimisation of code)
+"""
+import os
+import logging
+
+from eos.command_line import commandLineArgs
+
+def main():
+    logging.getLogger().setLevel(logging.CRITICAL)
+    clas = commandLineArgs([], 'eosls', extra_args=[
+        dict(dest='path', nargs='*', default=['.'], help='paths to list file in')])
+
+    from glob import glob
+    import tabulate
+    from eos.file_reader import AmorHeader
+
+    files = []
+    for path in clas.path:
+        files+=glob(os.path.join(path, 'amor*.hdf'))
+    files.sort()
+
+    data = {
+        'File name': [],
+        'Start Time': [],
+        'mu': [],
+        'nu': [],
+        'div': [],
+        'Sample': [],
+        'T [K]': [],
+        'H [T]': [],
+        }
+    for fi in files:
+        data['File name'].append(os.path.basename(fi))
+        ah = AmorHeader(fi)
+        data['Sample'].append(ah.sample.name)
+        data['Start Time'].append(ah.fileDate.strftime('%y %m-%d %H:%M:%S'))
+        data['mu'].append('%.3f' % ah.geometry.mu)
+        data['nu'].append('%.3f' % ah.geometry.nu)
+        data['div'].append('%.3f' % ah.geometry.div)
+
+        T = ah.sample.sample_parameters.get('temperature', None)
+        data['T [K]'].append(T.magnitude if T is not None else '-')
+
+        H = ah.sample.sample_parameters.get('magnetic_field', None)
+        data['H [T]'].append(H.magnitude if H is not None else '-')
+
+    print(tabulate.tabulate(data, headers="keys"))
+
+if __name__ == '__main__':
+    main()

eos/nicos.py

@@ -0,0 +1,46 @@
+"""
+amor-nicos vizualising data from Amor@SINQ, PSI
+
+Author: Jochen Stahn (algorithms, python draft),
+        Artur Glavic (structuring and optimisation of code)
+"""
+import logging
+
+# need to do absolute import here as pyinstaller requires it
+from eos.options import E2HConfig, ReaderConfig, ExperimentConfig, E2HReductionConfig
+from eos.command_line import commandLineArgs
+from eos.logconfig import setup_logging, update_loglevel
+
+
+def main():
+    setup_logging()
+    logging.getLogger('matplotlib').setLevel(logging.WARNING)
+
+    # read command line arguments and generate classes holding configuration parameters
+    clas = commandLineArgs([ReaderConfig, ExperimentConfig, E2HReductionConfig],
+                           'amor-nicos')
+    update_loglevel(clas.verbose)
+    if clas.verbose<2:
+        # only log info level in logfile
+        logger = logging.getLogger()  # logging.getLogger('quicknxs')
+        logger.setLevel(logging.INFO)
+
+    reader_config = ReaderConfig.from_args(clas)
+    experiment_config = ExperimentConfig.from_args(clas)
+    reduction_config = E2HReductionConfig.from_args(clas)
+    config = E2HConfig(reader_config, experiment_config, reduction_config)
+
+    logging.warning('######## amor-nicos - Nicos histogram for Amor ########')
+    from eos.reduction_kafka import KafkaReduction
+
+    # only import heavy module if sufficient command line parameters were provided
+    from eos.reduction_reflectivity import ReflectivityReduction
+    # Create reducer with these arguments
+    reducer = KafkaReduction(config)
+    # Perform actual reduction
+    reducer.reduce()
+
+    logging.info('######## amor-nicos - finished ########')
+
+if __name__ == '__main__':
+    main()

eos/nicos_interface.py

@@ -0,0 +1,60 @@
+"""
+Functions used to directly access information from nicos.
+"""
+
+import socket
+import platform
+import logging
+import subprocess
+
+ON_AMOR = platform.node().startswith('amor')
+NICOS_CACHE_DIR = '/home/data/nicosdata/cache/'
+GREP = '/usr/bin/grep "value"'
+
+
+def lookup_nicos_value(key, nicos_key, dtype=float, suffix='', year="2026"):
+    # TODO: Implement direct communication to nicos to read the cache
+    if nicos_key=='ignore':
+        return dtype(0)
+    try:
+        logging.debug(f'     trying socket request for device {nicos_key}')
+        response = nicos_single_request(nicos_key)
+        logging.info(f"     using parameter {nicos_key} from nicos cache via socket")
+        return dtype(response)
+    except Exception as e:
+        logging.debug(f'       socket request failed with {e!r}')
+    if ON_AMOR:
+        logging.debug(f"     trying to extract {nicos_key} from nicos cache files")
+        call = f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year}{suffix}'
+        try:
+            value = str(subprocess.getoutput(call)).split('\t')[-1]
+            logging.info(f"     using parameter {nicos_key} from nicos cache file")
+            return dtype(value)
+        except Exception:
+            logging.error(f"       couldn't get value from nicos cache {nicos_key}, {call}")
+            return dtype(0)
+    else:
+        logging.warning(f"     parameter {key} not found, relpace by zero")
+        return dtype(0)
+
+def nicos_single_request(device):
+    sentinel = b'\n'
+    with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s:
+        s.settimeout(1.0)
+        s.connect(('amor', 14869))
+
+        tosend = f'@nicos/{device}/value?\n'
+
+        # write request
+        # self.log.debug("get_explicit: sending %r", tosend)
+        s.sendall(tosend.encode())
+
+        # read response
+        data = b''
+        while not data.endswith(sentinel):
+            newdata = s.recv(8192)  # blocking read
+            if not newdata:
+                raise IOError('cache closed connection')
+            data += newdata
+        s.shutdown(socket.SHUT_RDWR)
+    return data.decode('utf-8').split('=')[-1]

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)
@@ -66,6 +67,30 @@ class LZNormalisation:
         self.monitor = 1.
         return self
 
+    @classmethod
+    def model(cls, grid:LZGrid) -> 'LZNormalisation':
+        # generate a normalization based on angular and wavelength distribution model
+        # TODO: add options for sample size for better absolute normalization
+        logging.warning(f'normalisation is model')
+        self = super().__new__(cls)
+        self.angle = 1.0
+        self.monitor = 4e6
+
+        lamda_l  = grid.lamda()
+        lamda_c = (lamda_l[:-1]+lamda_l[1:])/2
+
+        delta = np.rad2deg(np.arctan2(grid.z(), Detector.distance))/2.0
+        delta_c = (delta[:-1]+delta[1:])/2-delta.mean()
+
+        # approximate spectrum by Maxwell-Boltzmann and intensity by linear footprint
+        a = 3.8
+        Ilambda = np.sqrt(2./np.pi)*lamda_c**2/a**3*np.exp(-lamda_c**2/(2.*a**2))
+        Idelta = np.where(abs(delta_c)<0.75, (self.angle-delta_c), np.nan)
+
+        self.norm = 1e6*Ilambda[:, np.newaxis]*Idelta[np.newaxis, :]
+
+        return self
+
     def safe(self, filename, hash):
         with open(filename, 'wb') as fh:
             np.save(fh, hash, allow_pickle=False)
@@ -75,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

@@ -10,6 +10,7 @@ import numpy as np
 
 import logging
 
+
 try:
     from enum import StrEnum
 except ImportError:
@@ -20,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"
@@ -27,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):
         """
@@ -89,11 +96,14 @@ 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
+                    typ = get_args(typ)[0]
             elif get_origin(typ) is tuple:
                 args['nargs'] = len(get_args(typ))
                 typ = get_args(typ)[0]
-
             if issubclass(typ, StrEnum):
                 args['choices'] = [ci.value for ci in typ]
                 if field.default is not MISSING:
@@ -114,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
 
@@ -148,7 +159,12 @@ class ArgParsable:
             if get_origin(field.type) is Union and type(None) in get_args(field.type):
                 # optional argument
                 typ = get_args(field.type)[0]
-
+            if get_origin(typ) is list:
+                item_typ = get_args(typ)[0]
+                if get_origin(item_typ) is tuple:
+                    # tuple of items are put together during evaluation
+                    tuple_length = len(get_args(item_typ))
+                    value = [tuple(value[i*tuple_length+j] for j in range(tuple_length)) for i in range(len(value)//tuple_length)]
             if isinstance(typ, type) and issubclass(typ, StrEnum):
                 # convert str to enum
                 try:
@@ -160,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
@@ -170,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(
@@ -201,6 +246,14 @@ class MonitorType(StrEnum):
     neutron_monitor = 'n'
     debug = 'x'
 
+MONITOR_UNITS = {
+    MonitorType.neutron_monitor: 'cnts',
+    MonitorType.proton_charge: 'mC',
+    MonitorType.time: 's',
+    MonitorType.auto: 'various',
+    MonitorType.debug: 'mC',
+    }
+
 @dataclass
 class ExperimentConfig(ArgParsable):
     chopperPhase: float = field(
@@ -286,7 +339,7 @@ class ExperimentConfig(ArgParsable):
                 },                               
             )                          
     muOffset: Optional[float] = field(
-            default=0,
+            default=None,
             metadata={
                 'short': 'm',
                 'group': 'sample',  
@@ -308,7 +361,7 @@ class NormalisationMethod(StrEnum):
     under_illuminated = 'u'
 
 @dataclass
-class ReductionConfig(ArgParsable):
+class ReflectivityReductionConfig(ArgParsable):
     fileIdentifier: List[str] = field(
             metadata={
                 'short':    'f',
@@ -350,6 +403,14 @@ class ReductionConfig(ArgParsable):
                 'help': 'theta region of interest w.r.t. beam center',
                 },
             )
+    thetaFilters: List[Tuple[float, float]] = field(
+            default_factory=lambda: [],
+            metadata={
+                'short': 'TF',
+                'group': 'region of interest',
+                'help': 'add one or more theta ranges that will be filtered in reduction',
+                },
+            )
     normalisationMethod: NormalisationMethod = field(
             default=NormalisationMethod.over_illuminated,
             metadata={
@@ -357,6 +418,21 @@ class ReductionConfig(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={
@@ -380,7 +456,7 @@ class ReductionConfig(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,
@@ -395,32 +471,14 @@ class ReductionConfig(ArgParsable):
                 },
             )
 
-    def _expand_file_list(self, short_notation:str):
-        """Evaluate string entry for file number lists"""
-        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)]
-        file_list.sort()
-        return file_list
-
-    def data_files(self):
-        # get input files from expanding fileIdentifier
-        return list(map(self._expand_file_list, self.fileIdentifier))
-
-    def normal_files(self):
-        return list(map(self._expand_file_list, self.normalisationFileIdentifier))
+    logfilter: List[str] = field(
+            default_factory=lambda: [],
+            metadata={
+                'short': 'lf',
+                'group': 'region of interest',
+                'help':  'filter using comparison to a log values, multpiple allowd (example "sample_temperature<150")',
+                },
+            )
 
 
 class OutputFomatOption(StrEnum):
@@ -440,9 +498,10 @@ class PlotColormaps(StrEnum):
     inferno = "inferno"
     gist_rainbow = "gist_rainbow"
     nipy_spectral = "nipy_spectral"
+    jochen_deluxe = "jochen_deluxe"
 
 @dataclass
-class OutputConfig(ArgParsable):
+class ReflectivityOutputConfig(ArgParsable):
     outputFormats: List[OutputFomatOption] = field(
             default_factory=lambda: ['Rqz.ort'],
             metadata={
@@ -484,6 +543,14 @@ class OutputConfig(ArgParsable):
                 },
             )
 
+    append: bool = field(
+            default=False,
+            metadata={
+                'group': 'output',
+                'help': 'if file already exists, append result as additional ORSO dataset (only Rqz.ort)',
+                },
+            )
+
     def _output_format_list(self, outputFormat):
         format_list = []
         if OutputFomatOption.ort in outputFormat\
@@ -511,11 +578,11 @@ class OutputConfig(ArgParsable):
 # ===================================
 
 @dataclass
-class EOSConfig:
+class ReflectivityConfig:
     reader: ReaderConfig
     experiment: ExperimentConfig
-    reduction: ReductionConfig
-    output: OutputConfig
+    reduction: ReflectivityReductionConfig
+    output: ReflectivityOutputConfig
     
     _call_string_overwrite=None
     
@@ -528,84 +595,158 @@ class EOSConfig:
     
     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)}'
+        logging.debug(f'Argument list build in EOSConfig.call_string: {cpout}')
+        return  cpout
 
-        mlst = base + inpt + otpt 
-        if mask:
-            mlst += mask
-        if para:
-            mlst += para
-        if acts:
-            mlst += acts
-        if modl:
-            mlst += modl
+class E2HPlotSelection(StrEnum):
+    All = 'all'
+    Raw = 'raw'
+    YZ = 'Iyz'
+    LT = 'Ilt'
+    YT = 'Iyt'
+    TZ = 'Itz'
+    Q = 'Iq'
+    L = 'Il'
+    T = 'It'
+    ToF = 'tof'
 
-        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
+class E2HPlotArguments(StrEnum):
+    Default = 'def'
+    OutputFile = 'file'
+    Logarithmic = 'log'
+    Linear = 'lin'
 
-            
+@dataclass
+class E2HReductionConfig(ArgParsable):
+    fileIdentifier: str = field(
+            default='0',
+            metadata={
+                'short':    'f',
+                'priority': 100,
+                'group':    'input data',
+                'help':     'file number(s) or offset (if < 1), events2histogram only accepts one short code',
+                },
+            )
+
+    show_plot: bool = field(
+            default=False,
+            metadata={
+                'short': 'sp',
+                'group': 'output',
+                'help': 'show matplotlib graphs of results',
+                },
+            )
+
+    plot: E2HPlotSelection = field(
+            default=E2HPlotSelection.All,
+            metadata={
+                'short': 'p',
+                'group': 'output',
+                'help': 'select what to plot or write',
+                },
+            )
+
+    kafka: bool = field(
+            default=False,
+            metadata={
+                'group': 'output',
+                'help': 'send result to kafka for Nicos',
+                },
+            )
+
+    plotArgs: E2HPlotArguments = field(
+            default=E2HPlotArguments.Default,
+            metadata={
+                'short': 'pa',
+                'group': 'output',
+                'help': 'select configuration for plot',
+                },
+            )
+
+    update: bool = field(
+            default=False,
+            metadata={
+                'short': 'u',
+                'group': 'output',
+                'help': 'keep running in the background and update plot when file is modified, implies --plot',
+                },
+            )
+
+    fast: bool = field(
+            default=False,
+            metadata={
+                'group': 'input data',
+                'help': 'skip some reduction steps to speed up calculations',
+                },
+            )
+
+    normalizationModel: bool = field(
+            default=False,
+            metadata={
+                'short': 'nm',
+                'group': 'input data',
+                'help': 'use model for incoming spectrum and divergence to normalize for reflectivity',
+                },
+            )
+
+    plot_colormap: PlotColormaps = field(
+            default=PlotColormaps.jochen_deluxe,
+            metadata={
+                'short': 'pcmap',
+                'group': 'output',
+                'help': 'matplotlib colormap used in lambda-theta graphs when plotting',
+                },
+            )
+
+    max_events: int = field(
+            default = 10_000_000,
+            metadata={
+                'group': 'input data',
+                'help':  'maximum number of events read at once',
+                },
+            )
+
+    thetaRangeR: Tuple[float, float] = field(
+            default_factory=lambda: [-0.75, 0.75],
+            metadata={
+                'short': 'T',
+                'group': 'region of interest',
+                'help': 'theta region of interest w.r.t. beam center',
+                },
+            )
+
+    thetaFilters: List[Tuple[float, float]] = field(
+            default_factory=lambda: [],
+            metadata={
+                'short': 'TF',
+                'group': 'region of interest',
+                'help': 'add one or more theta ranges that will be filtered in reduction',
+                },
+            )
+
+    fontsize: float = field(
+            default=8.,
+            metadata={
+                'short': 'pf',
+                'group': 'output',
+                'help': 'font size for graphs',
+                },
+            )
+
+@dataclass
+class E2HConfig:
+    reader: ReaderConfig
+    experiment: ExperimentConfig
+    reduction: E2HReductionConfig

eos/path_handling.py

@@ -1,6 +1,7 @@
 """
 Defines how file paths are resolved from short_notation, year and number to filename.
 """
+import logging
 import os
 from typing import List
 
@@ -18,7 +19,7 @@ class PathResolver:
         """Evaluate string entry for file number lists"""
         file_list = []
         for i in short_notation.split(','):
-            if '-' in i:
+            if '-' in i and not i.startswith('-'):
                 if ':' in i:
                     step = i.split(':', 1)[1]
                     file_list += range(int(i.split('-', 1)[0]),
@@ -34,6 +35,8 @@ class PathResolver:
         return list(sorted(file_list))
 
     def get_path(self, number):
+        if number<=0:
+            number = self.search_latest(number)
         fileName = f'amor{self.year}n{number:06d}.hdf'
         path = ''
         for rawd in self.rawPath:
@@ -41,9 +44,39 @@ class PathResolver:
                 path = rawd
                 break
         if not path:
-            if os.path.exists(
-                    f'/afs/psi.ch/project/sinqdata/{self.year}/amor/{int(number/1000)}/{fileName}'):
-                path = f'/afs/psi.ch/project/sinqdata/{self.year}/amor/{int(number/1000)}'
+            from glob import glob
+            potential_file = glob(f'/home/amor/data/{self.year}/*/{fileName}')
+            if len(potential_file)>0:
+                path = os.path.dirname(potential_file[0])
             else:
-                raise FileNotFoundError(f'# ERROR: the file {fileName} can not be found in {self.rawPath}')
+                raise FileNotFoundError(f'# ERROR: the file {fileName} can not be found '
+                                        f'in {self.rawPath+["/home/amor/data"]}')
         return os.path.join(path, fileName)
+
+    def search_latest(self, number):
+        if number>0:
+            raise ValueError('number needs to be relative index (negative)')
+        if os.path.exists(f'/home/amor/data/{self.year}/DataNumber'):
+            try:
+                with open(f'/home/amor/data/{self.year}/DataNumber', 'r') as fh:
+                    current_index = int(fh.readline())-1
+            except Exception:
+                logging.error('Can not access DataNumber', exc_info=True)
+            else:
+                return current_index+number
+        # find all files from the given year, convert to number and return latest
+        from glob import glob
+        possible_files = []
+        for rawd in self.rawPath:
+            possible_files += glob(os.path.join(rawd, f'amor{self.year}n??????.hdf'))
+
+        possible_files += glob(f'/home/amor/data/{self.year}/*/amor{self.year}n??????.hdf')
+        possible_indices = list(set([int(os.path.basename(fi)[9:15]) for fi in possible_files]))
+        possible_indices.sort()
+        try:
+            return possible_indices[number-1]
+        except IndexError:
+            raise FileNotFoundError(f'# Could not find suitable file for relative index {number} '
+                                    f'in {self.rawPath+["/home/amor/data"]}, '
+                                    f'possible indices {possible_indices}')
+

eos/projection.py

@@ -3,6 +3,10 @@ Classes used to calculate projections/binnings from event data onto given grids.
 """
 
 import logging
+import warnings
+from abc import ABC, abstractmethod
+from typing import List, Tuple, Union
+
 import numpy as np
 from dataclasses import dataclass
 
@@ -10,6 +14,19 @@ from .event_data_types import EventDatasetProtocol
 from .instrument import Detector, LZGrid
 from .normalization import LZNormalisation
 
+class ProjectionInterface(ABC):
+    @abstractmethod
+    def project(self, dataset: EventDatasetProtocol, monitor: float): ...
+
+    @abstractmethod
+    def clear(self): ...
+
+    @abstractmethod
+    def plot(self, **kwargs): ...
+
+    @abstractmethod
+    def update_plot(self): ...
+
 @dataclass
 class ProjectedReflectivity:
     R: np.ndarray
@@ -70,22 +87,34 @@ class ProjectedReflectivity:
 
     def plot(self, **kwargs):
         from matplotlib import pyplot as plt
-        plt.errorbar(self.Q, self.R, xerr=self.dQ, yerr=self.dR, **kwargs)
+        plt.errorbar(self.Q, self.R, yerr=self.dR, **kwargs)
         plt.yscale('log')
-        plt.xlabel('Q / $\\AA^{-1}$')
+        plt.xlabel('Q / Å$^{-1}$')
         plt.ylabel('R')
 
-class LZProjection:
+class LZProjection(ProjectionInterface):
     grid: LZGrid
     lamda: np.ndarray
     alphaF: np.ndarray
     is_normalized: bool
+    angle: float
 
     data: np.recarray
+    _dtype = np.dtype([
+            ('I', np.float64),
+            ('mask', bool),
+            ('ref', np.float64),
+            ('err', np.float64),
+            ('res', np.float64),
+            ('qz', np.float64),
+            ('qx', np.float64),
+            ('norm', np.float64),
+            ])
 
     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()
@@ -95,16 +124,7 @@ class LZProjection:
 
         self.lamda  = lz_shape*lamda_c[:, np.newaxis]
         self.alphaF = lz_shape*alphaF_z[np.newaxis, :]
-        self.data = np.zeros(self.alphaF.shape, dtype=[
-            ('I', np.float64),
-            ('mask', bool),
-            ('ref', np.float64),
-            ('err', np.float64),
-            ('res', np.float64),
-            ('qz', np.float64),
-            ('qx', np.float64),
-            ('norm', np.float64),
-            ]).view(np.recarray)
+        self.data = np.zeros(self.alphaF.shape, dtype=self._dtype).view(np.recarray)
         self.data.mask = True
         self.monitor = 0.
 
@@ -156,15 +176,19 @@ class LZProjection:
         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):
         """
             Project dataset on grid and add to intensity.
             Can be called multiple times to sequentially add events.
         """
+        # TODO: maybe move monitor calculation in here instead of reduction?
         e = dataset.data.events
         int_lz, *_  = np.histogram2d(e.lamda, e.detZ, bins = (self.grid.lamda(), self.grid.z()))
         self.data.I += int_lz
@@ -172,6 +196,13 @@ class LZProjection:
         # in case the intensity changed one needs to normalize again
         self.is_normalized = False
 
+    def clear(self):
+        # empty data
+        self.data[:] = 0
+        self.data.mask = True
+        self.monitor = 0.
+        self.norm_monitor = 1.
+
     @property
     def I(self):
         output = self.data.I[:]
@@ -180,22 +211,26 @@ class LZProjection:
 
     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)
+        fp_corr_lz[fp_corr_lz<0] = 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
@@ -215,6 +250,7 @@ class LZProjection:
             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:
@@ -222,75 +258,27 @@ class LZProjection:
         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)
 
-    ############## potential speedup not used right now, needs to be tested ####################
-    @classmethod
-    def histogram2d_lz(cls, 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 = cls.devide_bin(lamda_e, detZ_e.astype(np.int64), bins[0], dimension)
-        return np.array(binning), bins[0], bins[1]
-
-    @classmethod
-    def devide_bin(cls, 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=np.int64).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 = cls.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 = cls.devide_bin(lambda_e[right_region], position_e[right_region],
-                                              lamda_edges[split_idx:], dimension)
-        return left_list+right_list
-
     def plot(self, **kwargs):
         from matplotlib import pyplot as plt
         from matplotlib.colors import LogNorm
@@ -302,16 +290,543 @@ class LZProjection:
             cmap=False
 
         if self.is_normalized:
-            if not 'norm' in kwargs:
-                kwargs['norm'] = LogNorm(2e-3, 2.0)
-            plt.pcolormesh(self.lamda, self.alphaF, self.data.ref, **kwargs)
-            if cmap:
-                plt.colorbar(label='R')
+            I = self.data.ref
         else:
-            if not 'norm' in kwargs:
-                kwargs['norm'] = LogNorm()
-            plt.pcolormesh(self.lamda, self.alphaF, self.data.I, **kwargs)
-            if cmap:
+            I = self.data.I
+
+
+        if not 'norm' in kwargs:
+            vmin = I[(I>0)].min()
+            vmax = np.nanmax(I)
+            kwargs['norm'] = LogNorm(vmin, vmax, clip=True)
+
+
+        # suppress warning for wrongly sorted y-axis pixels (blades overlap)
+        with warnings.catch_warnings(action='ignore', category=UserWarning):
+            self._graph = plt.pcolormesh(self.lamda, self.alphaF, I, **kwargs)
+        if cmap:
+            if self.is_normalized:
+                plt.colorbar(label='R')
+            else:
                 plt.colorbar(label='I / cpm')
         plt.xlabel('$\\lambda$ / $\\AA$')
         plt.ylabel('$\\Theta$ / °')
+        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')
+
+        self._graph_axis = plt.gca()
+        plt.connect('button_press_event', self.draw_qline)
+
+    def update_plot(self):
+        """
+        Inline update of previous plot by just updating the data.
+        """
+        from matplotlib.colors import LogNorm
+        if self.is_normalized:
+            I = self.data.ref
+        else:
+            I = self.data.I
+
+        if isinstance(self._graph.norm, LogNorm):
+            vmin = I[(I>0)].min()*0.5
+        else:
+            vmin = 0
+        vmax = np.nanmax(I)
+        self._graph.set_array(I)
+        self._graph.norm.vmin = vmin
+        self._graph.norm.vmax = vmax
+
+        if self.is_normalized:
+            self._graph.set_array(self.data.ref)
+        else:
+            self._graph.set_array(self.data.I)
+
+    def draw_qline(self, event):
+        if event.inaxes is not self._graph_axis:
+            return
+        from matplotlib import pyplot as plt
+        tbm = self._graph_axis.figure.canvas.manager.toolbar.mode
+        if event.button is plt.MouseButton.LEFT and tbm=='':
+            slope = event.ydata/event.xdata
+            xmax = 12.5
+            self._graph_axis.plot([0, xmax], [0, slope*xmax], '-', color='grey')
+            self._graph_axis.text(event.xdata, event.ydata, f'q={np.deg2rad(slope)*4.*np.pi:.3f}', backgroundcolor='white')
+            plt.draw()
+        if event.button is plt.MouseButton.RIGHT and tbm=='':
+            for art in list(self._graph_axis.lines)+list(self._graph_axis.texts):
+                art.remove()
+            plt.draw()
+
+ONLY_MAP = ['colorbar', 'cmap', 'norm']
+
+class ReflectivityProjector(ProjectionInterface):
+    lzprojection: LZProjection
+    data: ProjectedReflectivity
+    # TODO: maybe implement direct 1d projection in here
+
+    def __init__(self, lzprojection, norm):
+        self.lzprojection = lzprojection
+        self.norm = norm
+
+    def project(self, dataset: EventDatasetProtocol, monitor: float):
+        self.lzprojection.project(dataset, monitor)
+        self.lzprojection.normalize_over_illuminated(self.norm)
+        self.data = self.lzprojection.project_on_qz()
+
+    def clear(self):
+        self.lzprojection.clear()
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        for key in ONLY_MAP:
+            if key in kwargs: del(kwargs[key])
+
+        self._graph = plt.errorbar(self.data.Q, self.data.R, xerr=self.data.dQ, yerr=self.data.dR, **kwargs)
+        self._graph_axis = plt.gca()
+        plt.title('Reflectivity (might be improperly normalized)')
+        plt.yscale('log')
+        plt.xlabel('Q / $\\AA^{-1}$')
+        plt.ylabel('R')
+
+    def update_plot(self):
+        ln, _, (barsx, barsy) = self._graph
+
+        yerr_top = self.data.R+self.data.dR
+        yerr_bot = self.data.R-self.data.dR
+        xerr_top = self.data.Q+self.data.dQ
+        xerr_bot = self.data.Q-self.data.dQ
+
+        new_segments_x = [np.array([[xt, y], [xb, y]]) for xt, xb, y in zip(xerr_top, xerr_bot, self.data.R)]
+        new_segments_y = [np.array([[x, yt], [x, yb]]) for x, yt, yb in zip(self.data.Q, yerr_top, yerr_bot)]
+        barsx.set_segments(new_segments_x)
+        barsy.set_segments(new_segments_y)
+
+        ln.set_ydata(self.data.R)
+
+
+class YZProjection(ProjectionInterface):
+    y: np.ndarray
+    z: np.ndarray
+
+    data: np.recarray
+    _dtype = np.dtype([
+            ('cts', np.float64),
+            ('I', np.float64),
+            ('err', np.float64),
+            ])
+
+    def __init__(self):
+        self.z = np.arange(Detector.nBlades*Detector.nWires+1)-0.5
+        self.y = np.arange(Detector.nStripes+1)-0.5
+        self.data = np.zeros((self.y.shape[0]-1, self.z.shape[0]-1), dtype=self._dtype).view(np.recarray)
+        self.monitor = 0.
+
+    def project(self, dataset: EventDatasetProtocol, monitor: float):
+        detYi, detZi, detX, delta = Detector.pixelLookUp[dataset.data.events.pixelID-1].T
+
+        cts , *_ = np.histogram2d(detYi, detZi, bins=(self.y, self.z))
+        self.data.cts += cts
+        self.monitor += monitor
+
+        self.data.I = self.data.cts / self.monitor
+        self.data.err = np.sqrt(self.data.cts) / self.monitor
+
+    def clear(self):
+        self.data[:] = 0
+        self.monitor = 0.
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        from matplotlib.colors import LogNorm
+
+        if 'colorbar' in kwargs:
+            cmap=True
+            del(kwargs['colorbar'])
+        else:
+            cmap=False
+
+        vmax = self.data.I.max()
+
+        if not 'norm' in kwargs:
+            vmin = self.data.I[(self.data.I>0)].min()*0.5
+            kwargs['norm'] = LogNorm(vmin, vmax)
+
+        self._graph = plt.pcolormesh(self.y, self.z, self.data.I.T, **kwargs)
+        if cmap:
+            plt.colorbar(label='I / cpm')
+
+        plt.xlabel('Y')
+        plt.ylabel('Z')
+        plt.xlim(self.y[0], self.y[-1])
+        plt.ylim(self.z[-1], self.z[0])
+        plt.title('Horizontal Pixel vs. Vertical Pixel')
+
+        self._graph_axis = plt.gca()
+        plt.connect('button_press_event', self.draw_yzcross)
+
+    def update_plot(self):
+        """
+        Inline update of previous plot by just updating the data.
+        """
+        from matplotlib.colors import LogNorm
+        if isinstance(self._graph.norm, LogNorm):
+            vmin = self.data.I[(self.data.I>0)].min()*0.5
+        else:
+            vmin = 0
+        vmax = self.data.I.max()
+        self._graph.set_array(self.data.I.T)
+        self._graph.norm.vmin = vmin
+        self._graph.norm.vmax = vmax
+
+    def draw_yzcross(self, event):
+        if event.inaxes is not self._graph_axis:
+            return
+        from matplotlib import pyplot as plt
+        tbm = self._graph_axis.figure.canvas.manager.toolbar.mode
+        if event.button is plt.MouseButton.LEFT and tbm=='':
+            self._graph_axis.plot([event.xdata, event.xdata], [self.z[0], self.z[-1]], '-', color='grey')
+            self._graph_axis.plot([self.y[0], self.y[-1]], [event.ydata, event.ydata], '-', color='grey')
+            self._graph_axis.text(event.xdata, event.ydata, f'({event.xdata:.1f}, {event.ydata:.1f})', backgroundcolor='white')
+            plt.draw()
+        if event.button is plt.MouseButton.RIGHT and tbm=='':
+            for art in list(self._graph_axis.lines)+list(self._graph_axis.texts):
+                art.remove()
+            plt.draw()
+
+class YTProjection(YZProjection):
+    theta: np.ndarray
+
+    def __init__(self, tthh: float):
+        dd = Detector.delta_z[1]-Detector.delta_z[0]
+        delta = np.hstack([Detector.delta_z, Detector.delta_z[-1]+dd])-dd/2.
+        self.theta  = tthh + delta
+        super().__init__()
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        from matplotlib.colors import LogNorm
+
+        if 'colorbar' in kwargs:
+            cmap=True
+            del(kwargs['colorbar'])
+        else:
+            cmap=False
+
+        if not 'norm' in kwargs:
+            kwargs['norm'] = LogNorm()
+
+        self._graph = plt.pcolormesh(self.y, self.theta, self.data.I.T, **kwargs)
+        if cmap:
+            plt.colorbar(label='I / cpm')
+
+        plt.xlabel('Y')
+        plt.ylabel('Theta / °')
+        plt.xlim(self.y[0], self.y[-1])
+        plt.ylim(self.theta[-1], self.theta[0])
+        plt.title('Horizontal Pixel vs. Angle')
+
+        self._graph_axis = plt.gca()
+        plt.connect('button_press_event', self.draw_tzcross)
+
+    def draw_tzcross(self, event):
+        if event.inaxes is not self._graph_axis:
+            return
+        from matplotlib import pyplot as plt
+        tbm = self._graph_axis.figure.canvas.manager.toolbar.mode
+        if event.button is plt.MouseButton.LEFT and tbm=='':
+            self._graph_axis.plot([event.xdata, event.xdata], [self.theta[0], self.theta[-1]], '-', color='grey')
+            self._graph_axis.plot([self.y[0], self.y[-1]], [event.ydata, event.ydata], '-', color='grey')
+            self._graph_axis.text(event.xdata, event.ydata, f'({event.xdata:.1f}, {event.ydata:.1f})', backgroundcolor='white')
+            plt.draw()
+        if event.button is plt.MouseButton.RIGHT and tbm=='':
+            for art in list(self._graph_axis.lines)+list(self._graph_axis.texts):
+                art.remove()
+            plt.draw()
+
+
+class TofZProjection(ProjectionInterface):
+    tof: np.ndarray
+    z: np.ndarray
+
+    data: np.recarray
+    _dtype = np.dtype([
+            ('cts', np.float64),
+            ('I', np.float64),
+            ('err', np.float64),
+            ])
+
+    def __init__(self, tau, foldback=False, combine=1):
+        self.z = np.arange(Detector.nBlades*Detector.nWires+1)-0.5
+        if foldback:
+            self.tof = np.arange(0, tau, 0.0005*combine)
+        else:
+            self.tof = np.arange(0, 2*tau, 0.0005*combine)
+        self.data = np.zeros((self.tof.shape[0]-1, self.z.shape[0]-1), dtype=self._dtype).view(np.recarray)
+        self.monitor = 0.
+
+    def project(self, dataset: EventDatasetProtocol, monitor: float):
+        detYi, detZi, detX, delta = Detector.pixelLookUp[dataset.data.events.pixelID-1].T
+
+        cts , *_ = np.histogram2d(dataset.data.events.tof, detZi, bins=(self.tof, self.z))
+        self.data.cts += cts
+        self.monitor += monitor
+
+        self.data.I = self.data.cts / self.monitor
+        self.data.err = np.sqrt(self.data.cts) / self.monitor
+
+    def clear(self):
+        self.data[:] = 0
+        self.monitor = 0.
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        from matplotlib.colors import LogNorm
+
+        if 'colorbar' in kwargs:
+            cmap=True
+            del(kwargs['colorbar'])
+        else:
+            cmap=False
+
+        if not 'norm' in kwargs:
+            kwargs['norm'] = LogNorm()
+
+        self._graph = plt.pcolormesh(self.tof*1e3, self.z, self.data.I.T, **kwargs)
+        if cmap:
+            plt.colorbar(label='I / cpm')
+
+        plt.xlabel('Time of Flight / ms')
+        plt.ylabel('Z')
+        plt.xlim(self.tof[0]*1e3, self.tof[-1]*1e3)
+        plt.ylim(self.z[-1], self.z[0])
+        plt.title('Time of Flight vs. Vertical Pixel')
+
+        self._graph_axis = plt.gca()
+        plt.connect('button_press_event', self.draw_tzcross)
+
+    def update_plot(self):
+        """
+        Inline update of previous plot by just updating the data.
+        """
+        from matplotlib.colors import LogNorm
+        if isinstance(self._graph.norm, LogNorm):
+            vmin = self.data.I[(self.data.I>0)].min()*0.5
+        else:
+            vmin = 0
+        vmax = self.data.I.max()
+        self._graph.set_array(self.data.I.T)
+        self._graph.norm.vmin = vmin
+        self._graph.norm.vmax = vmax
+
+    def draw_tzcross(self, event):
+        if event.inaxes is not self._graph_axis:
+            return
+        from matplotlib import pyplot as plt
+        tbm = self._graph_axis.figure.canvas.manager.toolbar.mode
+        if event.button is plt.MouseButton.LEFT and tbm=='':
+            self._graph_axis.plot([event.xdata, event.xdata], [self.z[0], self.z[-1]], '-', color='grey')
+            self._graph_axis.plot([self.tof[0]*1e3, self.tof[-1]*1e3], [event.ydata, event.ydata], '-', color='grey')
+            self._graph_axis.text(event.xdata, event.ydata, f'({event.xdata:.2f}, {event.ydata:.1f})', backgroundcolor='white')
+            plt.draw()
+        if event.button is plt.MouseButton.RIGHT and tbm=='':
+            for art in list(self._graph_axis.lines)+list(self._graph_axis.texts):
+                art.remove()
+            plt.draw()
+
+class TofProjection(ProjectionInterface):
+    tof: np.ndarray
+
+    data: np.recarray
+    _dtype = np.dtype([
+            ('cts', np.float64),
+            ('I', np.float64),
+            ('err', np.float64),
+            ])
+
+    def __init__(self, tau, foldback=False):
+        if foldback:
+            self.tof = np.arange(0, tau, 0.0005)
+        else:
+            self.tof = np.arange(0, 2*tau, 0.0005)
+        self.data = np.zeros(self.tof.shape[0]-1, dtype=self._dtype).view(np.recarray)
+        self.monitor = 0.
+
+    def project(self, dataset: EventDatasetProtocol, monitor: float):
+        cts , *_ = np.histogram(dataset.data.events.tof, bins=self.tof)
+        self.data.cts += cts
+        self.monitor += monitor
+
+        self.data.I = self.data.cts / self.monitor
+        self.data.err = np.sqrt(self.data.cts) / self.monitor
+
+    def clear(self):
+        self.data[:] = 0
+        self.monitor = 0.
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        for key in ONLY_MAP:
+            if key in kwargs: del(kwargs[key])
+
+
+        self._graph = plt.plot(self.tof[:-1]*1e3, self.data.I, **kwargs)
+
+        plt.xlabel('Time of Flight / ms')
+        plt.ylabel('I / cpm')
+        plt.xlim(self.tof[0]*1e3, self.tof[-1]*1e3)
+        plt.title('Time of Flight')
+
+    def update_plot(self):
+        """
+        Inline update of previous plot by just updating the data.
+        """
+        self._graph[0].set_ydata(self.data.I.T)
+
+class LProjection(ProjectionInterface):
+    lamda: np.ndarray
+
+    data: np.recarray
+    _dtype = np.dtype([
+            ('cts', np.float64),
+            ('I', np.float64),
+            ('err', np.float64),
+            ])
+
+    def __init__(self):
+        self.lamda = np.linspace(3.0, 12.0, 91)
+        self.data = np.zeros(self.lamda.shape[0]-1, dtype=self._dtype).view(np.recarray)
+        self.monitor = 0.
+
+    def project(self, dataset: EventDatasetProtocol, monitor: float):
+        cts , *_ = np.histogram(dataset.data.events.lamda, bins=self.lamda)
+        self.data.cts += cts
+        self.monitor += monitor
+
+        self.data.I = self.data.cts / self.monitor
+        self.data.err = np.sqrt(self.data.cts) / self.monitor
+
+    def clear(self):
+        self.data[:] = 0
+        self.monitor = 0.
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        for key in ONLY_MAP:
+            if key in kwargs: del(kwargs[key])
+
+
+        self._graph = plt.plot(self.lamda[:-1], self.data.I, **kwargs)
+
+        plt.xlabel('Wavelength / Angstrom')
+        plt.ylabel('I / cpm')
+        plt.xlim(self.lamda[0], self.lamda[-1])
+        plt.title('Wavelength')
+
+    def update_plot(self):
+        """
+        Inline update of previous plot by just updating the data.
+        """
+        self._graph[0].set_ydata(self.data.I.T)
+
+class TProjection(ProjectionInterface):
+    theta: np.ndarray
+    z: np.ndarray
+
+    data: np.recarray
+    _dtype = np.dtype([
+            ('cts', np.float64),
+            ('I', np.float64),
+            ('err', np.float64),
+            ])
+
+    def __init__(self, tthh):
+        self.z = np.arange(Detector.nBlades*Detector.nWires+1)-0.5
+        dd = Detector.delta_z[1]-Detector.delta_z[0]
+        delta = np.hstack([Detector.delta_z, Detector.delta_z[-1]+dd])-dd/2.
+        self.theta = tthh+delta
+        self.data = np.zeros(self.theta.shape[0]-1, dtype=self._dtype).view(np.recarray)
+        self.monitor = 0.
+
+    def project(self, dataset: EventDatasetProtocol, monitor: float):
+        detYi, detZi, detX, delta = Detector.pixelLookUp[dataset.data.events.pixelID-1].T
+
+        cts , *_ = np.histogram(detZi, bins=self.z)
+        self.data.cts += cts
+        self.monitor += monitor
+
+        self.data.I = self.data.cts / self.monitor
+        self.data.err = np.sqrt(self.data.cts) / self.monitor
+
+    def clear(self):
+        self.data[:] = 0
+        self.monitor = 0.
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        for key in ONLY_MAP:
+            if key in kwargs: del(kwargs[key])
+
+
+        self._graph = plt.plot(self.theta[:-1], self.data.I, **kwargs)
+
+        plt.xlabel('Reflection Angle / °')
+        plt.ylabel('I / cpm')
+        plt.xlim(self.theta[-1], self.theta[0])
+        plt.title('Theta')
+
+    def update_plot(self):
+        """
+        Inline update of previous plot by just updating the data.
+        """
+        self._graph[0].set_ydata(self.data.I.T)
+
+
+class CombinedProjection(ProjectionInterface):
+    """
+    Allows to put multiple projections together to conveniently generate combined graphs.
+    """
+    projections: List[ProjectionInterface]
+    projection_placements: List[Union[Tuple[int, int], Tuple[int, int, int, int]]]
+    grid_size: Tuple[int, int]
+
+
+    def __init__(self, grid_rows, grid_cols, projections, projection_placements):
+        self.projections = projections
+        self.projection_placements = projection_placements
+        self.grid_size = grid_rows, grid_cols
+
+    def project(self, dataset: EventDatasetProtocol, monitor: float):
+        for pi in self.projections:
+            pi.project(dataset, monitor)
+
+    def clear(self):
+        for pi in self.projections:
+            pi.clear()
+
+    def plot(self, **kwargs):
+        from matplotlib import pyplot as plt
+        fig = plt.gcf()
+        axs = fig.add_gridspec(self.grid_size[0], self.grid_size[1])
+        # axs = fig.add_gridspec(self.grid_size[0]+1, self.grid_size[1],
+        #                        height_ratios=[1.0 for i in range(self.grid_size[0])]+[0.2])
+        self._axes = []
+        for pi, placement in zip(self.projections, self.projection_placements):
+            if len(placement) == 2:
+                ax = fig.add_subplot(axs[placement[0], placement[1]])
+            else:
+                ax = fig.add_subplot(axs[placement[0]:placement[1], placement[2]:placement[3]])
+            pi.plot(**dict(kwargs))
+        # Create the RangeSlider
+        # from matplotlib.widgets import RangeSlider
+        # slider_ax = fig.add_subplot(axs[self.grid_size[0], :])
+        # self._slider = RangeSlider(slider_ax, "Plot Range", 0., 1., valinit=(0., 1.))
+        # self._slider.on_changed(self.update_range)
+
+    def update_plot(self):
+        for pi in self.projections:
+            pi.update_plot()
+
+    # def update_range(self, event):
+    #     ...

eos/reduction_e2h.py

@@ -0,0 +1,318 @@
+"""
+Events 2 histogram, quick reduction of single file to display during experiment.
+Can be used as a live preview with automatic update when files are modified.
+"""
+
+import logging
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.colors import ListedColormap
+
+from orsopy import fileio
+from datetime import datetime
+
+from .file_reader import AmorEventData, AmorHeader
+from .header import Header
+from .instrument import LZGrid
+from .normalization import LZNormalisation
+from .options import E2HConfig, E2HPlotArguments, IncidentAngle, MonitorType, E2HPlotSelection
+from . import event_handling as eh
+from .path_handling import PathResolver
+from .projection import CombinedProjection, LProjection, LZProjection, ProjectionInterface, ReflectivityProjector, \
+    TofProjection,  TofZProjection,  TProjection, YTProjection, YZProjection
+
+NEEDS_LAMDA = (E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q, E2HPlotSelection.L)
+
+class E2HReduction:
+    config: E2HConfig
+    header: Header
+    event_actions: eh.EventDataAction
+
+    _last_mtime = 0.
+    projection: ProjectionInterface
+
+    def __init__(self, config: E2HConfig):
+        self.config = config
+
+        self.header = Header()
+
+        self.fig = plt.figure()
+        self.register_colormap()
+        self.prepare_actions()
+
+    def prepare_actions(self):
+        """
+        Does not do any actual reduction.
+        """
+        self.path_resolver = PathResolver(self.config.reader.year, self.config.reader.rawPath)
+        self.file_list = self.path_resolver.resolve(self.config.reduction.fileIdentifier)
+        self.file_index = 0
+        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
+
+        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
+
+        # 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:
+            logging.info('    Fast reduction always uses time normalization')
+            self.config.experiment.monitorType = MonitorType.time
+        self.event_actions |= eh.AssociatePulseWithMonitor(self.config.experiment.monitorType)
+        if self.config.experiment.monitorType in [MonitorType.proton_charge, MonitorType.debug]:
+            # the filtering only makes sense if using actual monitor data, not time
+            self.event_actions |= eh.FilterMonitorThreshold(self.config.experiment.lowCurrentThreshold)
+        if not self.config.reduction.fast:
+            self.event_actions |= eh.FilterStrangeTimes()
+            if self.config.reduction.plot in [E2HPlotSelection.YT, E2HPlotSelection.YZ]:
+                # perform time fold-back and apply yRange filter if not fast mode
+                self.event_actions |= ea.MergeFrames()
+                self.event_actions |= ea.AnalyzePixelIDs(self.config.experiment.yRange)
+            if self.config.reduction.plot==E2HPlotSelection.YT:
+                # 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 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)
+        self.event_actions |= eh.ApplyMask()
+
+        # plot dependant options
+        if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]:
+            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,
+                                          E2HPlotSelection.YZ, E2HPlotSelection.TZ]:
+            self.plot_kwds['colorbar'] = True
+            self.plot_kwds['cmap'] = str(self.config.reduction.plot_colormap)
+            if self.config.reduction.plotArgs==E2HPlotArguments.Linear:
+                self.plot_kwds['norm'] = None
+
+    def reduce(self):
+        if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]:
+            if self.config.reduction.normalizationModel:
+                self.norm = LZNormalisation.model(self.grid)
+            else:
+                self.norm = LZNormalisation.unity(self.grid)
+
+        self.prepare_graphs()
+
+        while self.file_index < len(self.file_list):
+            self.read_data()
+            self.add_data()
+
+        if self.config.reduction.plotArgs==E2HPlotArguments.OutputFile:
+            self.create_file_output()
+        if self.config.reduction.plotArgs!=E2HPlotArguments.OutputFile or self.config.reduction.show_plot:
+            self.create_graph()
+
+        if self.config.reduction.plotArgs==E2HPlotArguments.Default and not self.config.reduction.update:
+            # safe to image file if not auto-updating graph
+            plt.savefig(f'e2h_{self.config.reduction.plot}.png', dpi=300)
+        if self.config.reduction.kafka:
+            from .kafka_serializer import ESSSerializer
+            self.serializer = ESSSerializer()
+            self.fig.canvas.mpl_connect('close_event', self.serializer.end_command_thread)
+            self.serializer.start_command_thread()
+            self.serializer.send(self.projection)
+        if self.config.reduction.update:
+            self.timer = self.fig.canvas.new_timer(1000)
+            self.timer.add_callback(self.update)
+            self.timer.start()
+        if self.config.reduction.show_plot:
+            plt.show()
+
+
+    def register_colormap(self):
+        cmap = plt.colormaps['turbo'](np.arange(256))
+        cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1])
+        cmap = ListedColormap(cmap, name='jochen_deluxe', N=cmap.shape[0])
+        #cmap.set_bad((1.,1.,0.9))
+        plt.colormaps.register(cmap)
+
+    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'):
+            # adjust range based on detector center
+            thetaRange = [ti+tthh for ti in self.config.reduction.thetaRangeR]
+        else:
+            thetaRange = [tthh - last_file_header.geometry.div/2, tthh + last_file_header.geometry.div/2]
+
+        if self.config.reduction.plot==E2HPlotSelection.LT:
+            self.projection = LZProjection(tthh, self.grid)
+            if not self.config.reduction.fast:
+                self.projection.correct_gravity(last_file_header.geometry.detectorDistance)
+            self.projection.apply_lamda_mask(self.config.experiment.lambdaRange)
+            self.projection.apply_theta_mask(thetaRange)
+            for thi in self.config.reduction.thetaFilters:
+                self.projection.apply_theta_filter((thi[0]+tthh, thi[1]+tthh))
+            self.projection.apply_norm_mask(self.norm)
+
+        if self.config.reduction.plot==E2HPlotSelection.Q:
+            plz = LZProjection(tthh, self.grid)
+            if not self.config.reduction.fast:
+                plz.correct_gravity(last_file_header.geometry.detectorDistance)
+            plz.calculate_q()
+            plz.apply_lamda_mask(self.config.experiment.lambdaRange)
+            plz.apply_theta_mask(thetaRange)
+            for thi in self.config.reduction.thetaFilters:
+                self.projection.apply_theta_filter((thi[0]+tthh, thi[1]+tthh))
+            plz.apply_norm_mask(self.norm)
+            self.projection = ReflectivityProjector(plz, self.norm)
+
+        if self.config.reduction.plot==E2HPlotSelection.YZ:
+            self.projection = YZProjection()
+
+        if self.config.reduction.plot==E2HPlotSelection.YT:
+            self.projection = YTProjection(tthh)
+
+        if self.config.reduction.plot==E2HPlotSelection.T:
+            self.projection = TProjection(tthh)
+
+        if self.config.reduction.plot==E2HPlotSelection.L:
+            self.projection = LProjection()
+
+        if self.config.reduction.plot==E2HPlotSelection.TZ:
+            self.projection = TofZProjection(last_file_header.timing.tau, foldback=not self.config.reduction.fast)
+
+        if self.config.reduction.plot==E2HPlotSelection.ToF:
+            self.projection = TofProjection(last_file_header.timing.tau, foldback=not self.config.reduction.fast)
+
+        if self.config.reduction.plot==E2HPlotSelection.All:
+            plz = LZProjection(tthh, self.grid)
+            if not self.config.reduction.fast:
+                plz.correct_gravity(last_file_header.geometry.detectorDistance)
+            plz.calculate_q()
+            plz.apply_lamda_mask(self.config.experiment.lambdaRange)
+            plz.apply_theta_mask(thetaRange)
+            for thi in self.config.reduction.thetaFilters:
+                plz.apply_theta_filter((thi[0]+tthh, thi[1]+tthh))
+            plz.apply_norm_mask(self.norm)
+            pr = ReflectivityProjector(plz, self.norm)
+            pyz = YZProjection()
+            self.projection = CombinedProjection(3, 2, [plz, pyz, pr],
+                                                 [(0, 2, 0, 1), (0, 2, 1, 2), (2, 3, 0, 2)])
+
+        if self.config.reduction.plot==E2HPlotSelection.Raw:
+            del(self.plot_kwds['colorbar'])
+            # A combined graph that does not require longer calculations
+            plyt = YTProjection(tthh)
+            pltofz = TofZProjection(last_file_header.timing.tau, foldback=not self.config.reduction.fast)
+            pltof =  TofProjection(last_file_header.timing.tau, foldback=not self.config.reduction.fast)
+            plt = TProjection(tthh)
+
+            self.projection = CombinedProjection(3, 3, [plyt, pltofz, plt, pltof],
+                                                 [(0,2, 0, 1), (0, 2, 1, 3), (2,3, 0,1),(2,3,1,3)])
+
+    def read_data(self):
+        fileName = self.file_list[self.file_index]
+        self.dataset = AmorEventData(fileName, max_events=self.config.reduction.max_events)
+        if self.dataset.EOF or fileName==self.file_list[-1]:
+            self.file_index += 1
+        self.event_actions(self.dataset)
+        self.dataset.update_header(self.header)
+
+        self.header.measurement_data_files.append(fileio.File(file=os.path.basename(fileName),
+                                                              timestamp=self.dataset.fileDate))
+
+    def add_data(self):
+        self.monitor = self.dataset.data.pulses.monitor.sum()
+        self.projection.project(self.dataset, monitor=self.monitor)
+        if self.config.reduction.plot==E2HPlotSelection.LT:
+            self.projection.normalize_over_illuminated(self.norm)
+
+    def create_file_output(self):
+        raise NotImplementedError("Export to text output not yet implemented")
+
+    def create_title(self):
+        output = "Events to Histogram - "
+        output += ",".join(["#"+os.path.basename(fi)[9:15].lstrip('0') for fi in self.file_list])
+        output += f" ($\\mu$={self.dataset.geometry.mu:.2f} ;"
+        output += f" $\\nu$={self.dataset.geometry.nu:.2f})"
+        if self.config.reduction.update:
+            output += f"\n at "+datetime.now().strftime("%m/%d/%Y %H:%M:%S")
+        return output
+
+    def create_graph(self):
+        plt.suptitle(self.create_title())
+        self.projection.plot(**self.plot_kwds)
+        plt.tight_layout(pad=0.5)
+
+    def replace_dataset(self, latest):
+        new_files = self.path_resolver.resolve(f'{latest}')
+        if not os.path.exists(new_files[-1]):
+            return
+        try:
+            # check that events exist in the new file
+            AmorEventData(new_files[-1], 0, max_events=1_000)
+        except Exception:
+            logging.debug("Problem when trying to load new dataset", exc_info=True)
+            return
+
+        logging.warning(f"Preceding to next file {latest}")
+        self.file_list = new_files
+        self.file_index = 0
+        self.prepare_actions()
+        self.prepare_graphs()
+        self.read_data()
+        self.add_data()
+        self.fig.clear()
+        self.create_graph()
+        plt.draw()
+
+    def update(self):
+        logging.debug("    check for update")
+        if self.config.reduction.fileIdentifier=='0':
+            # if latest file was choosen, check if new one available and switch to it
+            current = int(os.path.basename(self.file_list[-1])[9:15])
+            latest = self.path_resolver.search_latest(0)
+            if latest>current:
+                self.replace_dataset(latest)
+                return
+        # if all events were read last time, only load more if file was modified
+        if self.dataset.EOF and os.path.getmtime(self.file_list[-1])<=self._last_mtime:
+            return
+
+        self._last_mtime = os.path.getmtime(self.file_list[-1])
+        try:
+            update_data = AmorEventData(self.file_list[-1], self.dataset.last_index+1,
+                                        max_events=self.config.reduction.max_events)
+        except EOFError:
+            return
+        logging.info("    updating with new data")
+
+        self.event_actions(update_data)
+        self.dataset=update_data
+        self.monitor = self.dataset.data.pulses.monitor.sum()
+        self.projection.project(update_data, self.monitor)
+
+        self.projection.update_plot()
+        plt.suptitle(self.create_title())
+        plt.draw()
+
+        if self.config.reduction.kafka:
+            self.serializer.send(self.projection)

eos/reduction_kafka.py

@@ -0,0 +1,128 @@
+"""
+Events 2 histogram, quick reduction of single file to display during experiment.
+Can be used as a live preview with automatic update when files are modified.
+"""
+
+import logging
+import os
+
+from time import sleep
+from .kafka_events import KafkaEventData
+from .header import Header
+from .options import E2HConfig
+from . import event_handling as eh, event_analysis as ea
+from .projection import  TofZProjection,  YZProjection
+from .kafka_serializer import ESSSerializer
+
+
+class KafkaReduction:
+    config: E2HConfig
+    header: Header
+    event_actions: eh.EventDataAction
+
+    _last_mtime = 0.
+    proj_yz: YZProjection
+    proj_tofz = TofZProjection
+
+    def __init__(self, config: E2HConfig):
+        self.config = config
+
+        self.header = Header()
+        self.event_data = KafkaEventData()
+        self.event_data.start()
+
+        self.prepare_actions()
+
+    def prepare_actions(self):
+        """
+        Does not do any actual reduction.
+        """
+        # Actions on datasets not used for normalization
+        self.event_actions = eh.ApplyPhaseOffset(self.config.experiment.chopperPhaseOffset)
+        self.event_actions |= eh.CorrectChopperPhase()
+        self.event_actions |= ea.MergeFrames()
+        self.event_actions |= eh.ApplyMask()
+
+    def reduce(self):
+        self.create_projections()
+        self.read_data()
+        self.add_data()
+
+        self.serializer = ESSSerializer()
+        self.serializer.start_command_thread()
+
+        self.loop()
+
+    def create_projections(self):
+        self.proj_yz = YZProjection()
+        self.proj_tofz = TofZProjection(self.event_data.timing.tau, foldback=True, combine=2)
+
+    def read_data(self):
+        # make sure the first events have arrived before starting analysis
+        self.event_data.new_events.wait()
+        self.dataset = self.event_data.get_events()
+        self.event_actions(self.dataset)
+
+
+    def add_data(self):
+        self.monitor = self.dataset.monitor
+        self.proj_yz.project(self.dataset, monitor=self.monitor)
+        self.proj_tofz.project(self.dataset, monitor=self.monitor)
+
+    def loop(self):
+        self.wait_for = self.serializer.new_count_started
+        while True:
+            try:
+                self.update()
+                self.wait_for.wait(1.0)
+            except KeyboardInterrupt:
+                self.event_data.stop_event.set()
+                self.event_data.join()
+                self.serializer.end_command_thread()
+                return
+
+    def update(self):
+        if self.serializer.new_count_started.is_set():
+            logging.warning('Start new count, clearing event data')
+            self.wait_for = self.serializer.count_stopped
+            self.event_data.restart()
+            self.serializer.new_count_started.clear()
+            self.create_projections()
+            return
+        elif self.serializer.count_stopped.is_set() and not self.event_data.stop_counting.is_set():
+            return self.finish_count()
+        try:
+            update_data = self.event_data.get_events()
+        except EOFError:
+            return
+        logging.info("    updating with new data")
+
+        self.event_actions(update_data)
+        self.dataset=update_data
+        self.monitor = self.dataset.monitor
+        self.proj_yz.project(update_data, self.monitor)
+        self.proj_tofz.project(update_data, self.monitor)
+
+        self.serializer.send(self.proj_yz)
+        self.serializer.send(self.proj_tofz)
+
+    def finish_count(self):
+        logging.debug("    stop event set, hold event collection and send final results")
+        self.wait_for = self.serializer.new_count_started
+        self.event_data.stop_counting.set()
+
+        try:
+            update_data = self.event_data.get_events()
+        except EOFError:
+            pass
+        else:
+            self.event_actions(update_data)
+            self.dataset = update_data
+            self.monitor = self.dataset.monitor
+            self.proj_yz.project(update_data, self.monitor)
+            self.proj_tofz.project(update_data, self.monitor)
+
+        logging.warning(f'  stop counting, total events {int(self.proj_tofz.data.cts.sum())}')
+
+        self.serializer.send(self.proj_yz, final=True)
+        self.serializer.send(self.proj_tofz, final=True)

eos/reduction_reflectivity.py

@@ -5,30 +5,26 @@ import sys
 import numpy as np
 from orsopy import fileio
 
+from .event_analysis import FilterByLog
+from .event_handling import ApplyMask
 from .file_reader import AmorEventData
 from .header import Header
 from .path_handling import PathResolver
-from .options import EOSConfig, IncidentAngle, MonitorType, NormalisationMethod
-from .instrument import Detector, LZGrid
+from .options import ReflectivityConfig, IncidentAngle, MonitorType, NormalisationMethod, MONITOR_UNITS
+from .instrument import LZGrid
 from .normalization import LZNormalisation
 from . import event_handling as eh, event_analysis as ea
 from .projection import LZProjection
 
 
-MONITOR_UNITS = {
-    MonitorType.neutron_monitor: 'cnts',
-    MonitorType.proton_charge: 'mC',
-    MonitorType.time: 's',
-    MonitorType.auto: 'various',
-    MonitorType.debug: 'mC',
-    }
-
-class AmorReduction:
-    def __init__(self, config: EOSConfig):
-        self.experiment_config = config.experiment
-        self.reader_config = config.reader
-        self.reduction_config = config.reduction
-        self.output_config = config.output
+class ReflectivityReduction:
+    config: ReflectivityConfig
+    header: Header
+    normevent_actions: eh.EventDataAction
+    dataevent_actions: eh.EventDataAction
+    
+    def __init__(self, config: ReflectivityConfig):
+        self.config = config
 
         self.header = Header()
         self.header.reduction.call = config.call_string()
@@ -37,61 +33,69 @@ class AmorReduction:
 
     def prepare_actions(self):
         """
-        Does not do any actual reduction.
+        Prepare the actions applied to each event dataset, does not do any actual reduction.
         """
-        self.path_resolver = PathResolver(self.reader_config.year, self.reader_config.rawPath)
+        self.path_resolver = PathResolver(self.config.reader.year, self.config.reader.rawPath)
 
         # setup all actions performed on event datasets before projection on the grid
         # The order of these corrections matter as some rely on parameters modified before
-        if self.reduction_config.normalisationFileIdentifier:
+        if self.config.reduction.normalisationFileIdentifier:
             # explicit steps performed on AmorEventDataset for normalization files
-            self.normevent_actions = eh.ApplyPhaseOffset(self.experiment_config.chopperPhaseOffset)
+            self.normevent_actions = eh.ApplyPhaseOffset(self.config.experiment.chopperPhaseOffset)
             self.normevent_actions |= eh.CorrectChopperPhase()
-            if self.experiment_config.monitorType in [MonitorType.proton_charge, MonitorType.debug]:
+            self.normevent_actions |= eh.AssociatePulseWithMonitor(self.config.experiment.monitorType)
+            if self.config.experiment.monitorType in [MonitorType.proton_charge, MonitorType.debug]:
                 self.normevent_actions |= ea.ExtractWalltime()
-            self.normevent_actions |= eh.AssociatePulseWithMonitor(self.experiment_config.monitorType,
-                                                                   self.experiment_config.lowCurrentThreshold)
+                self.normevent_actions |= eh.FilterMonitorThreshold(self.config.experiment.lowCurrentThreshold)
             self.normevent_actions |= eh.FilterStrangeTimes()
             self.normevent_actions |= ea.MergeFrames()
-            self.normevent_actions |= ea.AnalyzePixelIDs(self.experiment_config.yRange)
-            self.normevent_actions |= eh.TofTimeCorrection(self.experiment_config.incidentAngle==IncidentAngle.alphaF)
-            self.normevent_actions |= ea.CalculateWavelength(self.experiment_config.lambdaRange)
+            self.normevent_actions |= ea.AnalyzePixelIDs(self.config.experiment.yRange)
+            self.normevent_actions |= eh.TofTimeCorrection(self.config.experiment.incidentAngle==IncidentAngle.alphaF)
+            self.normevent_actions |= ea.CalculateWavelength(self.config.experiment.lambdaRange)
             self.normevent_actions |= eh.ApplyMask()
         # Actions on datasets not used for normalization
-        self.dataevent_actions = eh.ApplyPhaseOffset(self.experiment_config.chopperPhaseOffset)
-        self.dataevent_actions |= eh.ApplyParameterOverwrites(self.experiment_config) # some actions use instrument parameters, change before that
+        self.dataevent_actions = eh.ApplyPhaseOffset(self.config.experiment.chopperPhaseOffset)
+        self.dataevent_actions |= eh.ApplyParameterOverwrites(self.config.experiment) # some actions use instrument parameters, change before that
         self.dataevent_actions |= eh.CorrectChopperPhase()
         self.dataevent_actions |= ea.ExtractWalltime()
         self.dataevent_time_correction = eh.CorrectSeriesTime(0) # will be set from first dataset
         self.dataevent_actions |= self.dataevent_time_correction
-        self.dataevent_actions |= eh.AssociatePulseWithMonitor(self.experiment_config.monitorType,
-                                                               self.experiment_config.lowCurrentThreshold)
+        self.dataevent_actions |= eh.AssociatePulseWithMonitor(self.config.experiment.monitorType)
+        if self.config.experiment.monitorType in [MonitorType.proton_charge or MonitorType.debug]:
+            # the filtering only makes sense if using actual monitor data, not time
+            self.dataevent_actions |= eh.FilterMonitorThreshold(self.config.experiment.lowCurrentThreshold)
         self.dataevent_actions |= eh.FilterStrangeTimes()
         self.dataevent_actions |= ea.MergeFrames()
-        self.dataevent_actions |= ea.AnalyzePixelIDs(self.experiment_config.yRange)
-        self.dataevent_actions |= eh.TofTimeCorrection(self.experiment_config.incidentAngle==IncidentAngle.alphaF)
-        self.dataevent_actions |= ea.CalculateWavelength(self.experiment_config.lambdaRange)
-        self.dataevent_actions |= ea.CalculateQ(self.experiment_config.incidentAngle)
-        self.dataevent_actions |= ea.FilterQzRange(self.reduction_config.qzRange)
+        self.dataevent_actions |= ea.AnalyzePixelIDs(self.config.experiment.yRange)
+        self.dataevent_actions |= eh.TofTimeCorrection(self.config.experiment.incidentAngle==IncidentAngle.alphaF)
+        self.dataevent_actions |= ea.CalculateWavelength(self.config.experiment.lambdaRange)
+        self.dataevent_actions |= ea.CalculateQ(self.config.experiment.incidentAngle)
+        if not self.config.reduction.is_default('qzRange'):
+            self.dataevent_actions |= ea.FilterQzRange(self.config.reduction.qzRange)
+        for lf in self.config.reduction.logfilter:
+            self.dataevent_actions |= ea.FilterByLog(lf)
         self.dataevent_actions |= eh.ApplyMask()
 
-        self.grid = LZGrid(self.reduction_config.qResolution, self.reduction_config.qzRange)
+        self.grid = LZGrid(self.config.reduction.qResolution, self.config.reduction.qzRange)
 
     def reduce(self):
-        if not os.path.exists(f'{self.output_config.outputPath}'):
-            logging.debug(f'Creating destination path {self.output_config.outputPath}')
-            os.system(f'mkdir {self.output_config.outputPath}')
+        if not os.path.exists(f'{self.config.output.outputPath}'):
+            logging.debug(f'Creating destination path {self.config.output.outputPath}')
+            os.system(f'mkdir {self.config.output.outputPath}')
 
         # load or create normalisation matrix
-        if self.reduction_config.normalisationFileIdentifier:
+        if self.config.reduction.normalisationFileIdentifier:
             # TODO: change option definition to single normalization short_code
-            self.create_normalisation_map(self.reduction_config.normalisationFileIdentifier[0])
+            self.create_normalisation_map(self.config.reduction.normalisationFileIdentifier[0])
         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.reduction_config.subtract:
-            self.sq_q, self.sR_q, self.sdR_q, self.sFileName = self.loadRqz(self.reduction_config.subtract)
+        if self.config.reduction.subtract:
+            self.sq_q, self.sR_q, self.sdR_q, self.sFileName = self.loadRqz(self.config.reduction.subtract)
             logging.warning(f'loaded background file: {self.sFileName}')
             self.header.reduction.corrections.append(f'background from \'{self.sFileName}\' subtracted')
             self.subtract = True
@@ -101,44 +105,54 @@ class AmorReduction:
         # load measurement data and do the reduction
         self.datasetsRqz = []
         self.datasetsRlt = []
-        for i, short_notation in enumerate(self.reduction_config.fileIdentifier):
+        for i, short_notation in enumerate(self.config.reduction.fileIdentifier):
             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().replace(':', '-')
+            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.output_config.outputFormats:
+        if 'Rqz.ort' in self.config.output.outputFormats:
             self.save_Rqz()
 
-        if 'Rlt.ort' in self.output_config.outputFormats:
+        if 'Rlt.ort' in self.config.output.outputFormats:
             self.save_Rtl()
 
-        if self.output_config.plot:
+        if self.config.output.plot:
             import matplotlib.pyplot as plt
-            if 'Rqz.ort' in self.output_config.outputFormats:
+            if 'Rqz.ort' in self.config.output.outputFormats:
                 plt.figure(num=99)
                 plt.legend()
             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 = []
 
         self.dataset = AmorEventData(file_list[0])
-        if self.experiment_config.monitorType==MonitorType.auto:
+        if self.config.experiment.monitorType==MonitorType.auto:
             if self.dataset.data.proton_current.current.sum()>1:
-                self.experiment_config.monitorType = MonitorType.proton_charge
+                self.config.experiment.monitorType = MonitorType.proton_charge
                 logging.debug('      monitor type set to "proton current"')
             else:
-                self.experiment_config.monitorType = MonitorType.time
+                self.config.experiment.monitorType = MonitorType.time
                 logging.debug('      monitor type set to "time"')
             # update actions to sue selected monitor
             self.prepare_actions()
             # reload normalization to make sure the monitor matches
-            if self.reduction_config.normalisationFileIdentifier:
-                self.create_normalisation_map(self.reduction_config.normalisationFileIdentifier[0])
+            if self.config.reduction.normalisationFileIdentifier:
+                self.create_normalisation_map(self.config.reduction.normalisationFileIdentifier[0])
 
         self.dataevent_time_correction.seriesStartTime = self.dataset.eventStartTime
         self.dataevent_actions(self.dataset)
@@ -150,38 +164,62 @@ class AmorReduction:
             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))
 
-
-        if self.reduction_config.timeSlize:
+        if 'polarization_config_label' in self.dataset.data.device_logs:
+            pols = np.unique(self.dataset.data.device_logs['polarization_config_label'].value)
+            pols = pols[pols>0]
+            if len(pols)>1:
+                logging.warning(f'    found {len(pols)} polarization configurations, splitting dataset accordingly')
+                from copy import deepcopy
+                from . import const
+                full_ds = deepcopy(self.dataset)
+                for pi in pols:
+                    plabel = const.polarizationLabels[pi]
+                    pol_filter = FilterByLog(f'polarization_config_label=={pi}',
+                                             remove_switchpulse=True) | ApplyMask()
+                    logging.info(f'    filter {plabel} using polarization_config_label=={pi}')
+                    pol_filter(self.dataset)
+                    self.dataset.update_header(self.header)
+                    pol_filter.update_header(self.header)
+                    if self.config.reduction.timeSlize:
+                        if i>0:
+                            logging.warning(
+                                "    time slizing should only be used for one set of datafiles, check parameters")
+                        self.analyze_timeslices(i, polstr=f' : polarization = {plabel}')
+                    else:
+                        self.analyze_unsliced(i, polstr=f' : polarization = {plabel}')
+                    self.dataset = deepcopy(full_ds)
+                return
+        if self.config.reduction.timeSlize:
             if i>0:
                 logging.warning("    time slizing should only be used for one set of datafiles, check parameters")
             self.analyze_timeslices(i)
         else:
             self.analyze_unsliced(i)
 
-    def analyze_unsliced(self, i):
-        self.monitor = np.sum(self.dataset.data.pulses.monitor)
-        logging.warning(f'    monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.experiment_config.monitorType]}')
+    def analyze_unsliced(self, i, polstr=''):
+        self.monitor = self.dataset.data.pulses.monitor.sum()
+        logging.info(f'    monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.config.experiment.monitorType]}')
 
         proj:LZProjection = self.project_on_lz()
         try:
-            scale = self.reduction_config.scale[i]
+            scale = self.config.reduction.scale[i]
         except IndexError:
-            scale = self.reduction_config.scale[-1]
+            scale = self.config.reduction.scale[-1]
         proj.scale(scale)
 
-        if 'Rqz.ort' in self.output_config.outputFormats:
+        if 'Rqz.ort' in self.config.output.outputFormats:
             headerRqz = self.header.orso_header()
-            headerRqz.data_set = f'Nr {i} : mu = {self.dataset.geometry.mu:6.3f} deg'
+            headerRqz.data_set = f'Nr {i} : mu = {self.dataset.geometry.mu:6.3f} deg{polstr}'
 
             # projection on qz-grid
             result = proj.project_on_qz()
 
-            if self.reduction_config.autoscale:
+            if self.config.reduction.autoscale:
                 if i==0:
-                    result.autoscale(self.reduction_config.autoscale)
+                    result.autoscale(self.config.reduction.autoscale)
                 else:
                     result.stitch(self.last_result)
 
@@ -196,12 +234,12 @@ class AmorReduction:
             self.last_result = result
             self.datasetsRqz.append(orso_data)
 
-            if self.output_config.plot:
+            if self.config.output.plot:
                 import matplotlib.pyplot as plt
                 # plot all reflectivity results in same graph
                 plt.figure(num=99)
-                result.plot(label=f'{self.reduction_config.fileIdentifier[i]}')
-        if 'Rlt.ort' in self.output_config.outputFormats:
+                result.plot(label=f'{self.config.reduction.fileIdentifier[i]}')
+        if 'Rlt.ort' in self.config.output.outputFormats:
             columns = [
                 fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'),
                 fileio.Column('R', '', 'specular reflectivity'),
@@ -243,22 +281,22 @@ class AmorReduction:
                 self.datasetsRlt.append(orso_data)
                 j += 1
 
-            if self.output_config.plot:
+            if self.config.output.plot:
                 import matplotlib.pyplot as plt
                 plt.figure()
-                proj.plot(colorbar=True, cmap=str(self.output_config.plot_colormap))
-                plt.title(f'{self.reduction_config.fileIdentifier[i]}')
+                proj.plot(colorbar=True, cmap=str(self.config.output.plot_colormap))
+                plt.title(f'{self.config.reduction.fileIdentifier[i]}')
 
-    def analyze_timeslices(self, i):
+    def analyze_timeslices(self, i, polstr=''):
         wallTime_e = np.float64(self.dataset.data.events.wallTime)/1e9
         pulseTimeS = np.float64(self.dataset.data.pulses.time)/1e9
-        interval = self.reduction_config.timeSlize[0]
+        interval = self.config.reduction.timeSlize[0]
         try:
-            start = self.reduction_config.timeSlize[1]
+            start = self.config.reduction.timeSlize[1]
         except IndexError:
             start = 0
         try:
-            stop = self.reduction_config.timeSlize[2]
+            stop = self.config.reduction.timeSlize[2]
         except IndexError:
             stop = wallTime_e[-1]
         # make overwriting log lines possible by removing newline at the end
@@ -268,23 +306,23 @@ class AmorReduction:
         for ti, time in enumerate(np.arange(start, stop, interval)):
             slice = self.dataset.get_timeslice(time, time+interval)
             self.monitor = np.sum(slice.data.pulses.monitor)
-            logging.info(f'      {ti:<4d}  {time:6.0f}  {self.monitor:7.2f} {MONITOR_UNITS[self.experiment_config.monitorType]}')
+            logging.info(f'      {ti:<4d}  {time:6.0f}  {self.monitor:7.2f} {MONITOR_UNITS[self.config.experiment.monitorType]}')
 
             proj: LZProjection = self.project_on_lz(slice)
             try:
-                scale = self.reduction_config.scale[i]
+                scale = self.config.reduction.scale[i]
             except IndexError:
-                scale = self.reduction_config.scale[-1]
+                scale = self.config.reduction.scale[-1]
             proj.scale(scale)
 
             # projection on qz-grid
             result = proj.project_on_qz()
 
-            if self.reduction_config.autoscale:
+            if self.config.reduction.autoscale:
                 # scale every slice the same
                 if ti==0:
                     if i==0:
-                        atscale = result.autoscale(self.reduction_config.autoscale)
+                        atscale = result.autoscale(self.config.reduction.autoscale)
                     else:
                         atscale = result.stitch(self.last_result)
                 else:
@@ -299,15 +337,15 @@ class AmorReduction:
 
             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'
+            headerRqz.data_set = f'{i}_{ti}: time = {time:8.1f} s  to {time+interval:8.1f} s{polstr}'
             orso_data = fileio.OrsoDataset(headerRqz, result.data_for_time(time))
             self.datasetsRqz.append(orso_data)
 
-            if self.output_config.plot:
+            if self.config.output.plot:
                 import matplotlib.pyplot as plt
                 # plot all reflectivity results in same graph
                 plt.figure(num=99)
-                result.plot(label=f'{self.reduction_config.fileIdentifier[i]} @ {time:.1f}s')
+                result.plot(label=f'{self.config.reduction.fileIdentifier[i]} @ {time:.1f}s')
 
         self.last_result = result
         # reset normal logging behavior
@@ -315,19 +353,34 @@ class AmorReduction:
         logging.info(f'      done  {min(time+interval, pulseTimeS[-1]):5.0f}')
 
     def save_Rqz(self):
-        fname = os.path.join(self.output_config.outputPath, f'{self.output_config.outputName}.Rqz.ort')
+        fname = os.path.join(self.config.output.outputPath, f'{self.config.output.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)
+        if os.path.exists(fname) and self.config.output.append:
+            logging.info('        file already exists, append as new dataset')
+            with open(fname, 'r') as f:
+                f.readline()
+                theSecondLine = f.readline()[3:]
+            prev_data = fileio.load_orso(fname)
+            prev_names = [di.info.data_set for di in prev_data]
+            for i, di in enumerate(self.datasetsRqz):
+                while di.info.data_set in prev_names:
+                    if di.info.data_set.startswith('Nr '):
+                        di.info.data_set = f'Nr {i+len(prev_data)} :'+di.info.data_set.split(':', 1)[1]
+                        break
+                    di.info.data_set = di.info.data_set+'_'
+            fileio.save_orso(prev_data+self.datasetsRqz, fname, data_separator='\n', comment=theSecondLine)
+        else:
+            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.output_config.outputPath, f'{self.output_config.outputName}.Rlt.ort')
+        fname = os.path.join(self.config.output.outputPath, f'{self.config.output.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 loadRqz(self, name):
-        fname = os.path.join(self.output_config.outputPath, name)
+        fname = os.path.join(self.config.output.outputPath, name)
         if os.path.exists(fname):
             fileName = fname
         elif os.path.exists(f'{fname}.Rqz.ort'):
@@ -340,7 +393,7 @@ class AmorReduction:
         return q_q, Sq_q, dS_q, fileName
 
     def create_normalisation_map(self, short_notation):
-        outputPath = self.output_config.outputPath
+        outputPath = self.config.output.outputPath
         normalisation_list = self.path_resolver.expand_file_list(short_notation)
         name = '_'.join(map(str, normalisation_list))
         n_path = os.path.join(outputPath, f'{name}.norm')
@@ -364,44 +417,48 @@ class AmorReduction:
                 toadd = AmorEventData(nfi)
                 self.normevent_actions(toadd)
                 reference.append(toadd)
-            self.norm = LZNormalisation(reference, self.reduction_config.normalisationMethod, self.grid)
+            self.norm = LZNormalisation(reference, self.config.reduction.normalisationMethod, self.grid)
             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):
         if dataset is None:
             dataset=self.dataset
         proj = LZProjection.from_dataset(dataset, self.grid,
-                                         has_offspecular=(self.experiment_config.incidentAngle!=IncidentAngle.alphaF))
+                                         has_offspecular=(self.config.experiment.incidentAngle!=IncidentAngle.alphaF))
 
-        if not self.reduction_config.is_default('thetaRangeR'):
-            t0 = dataset.geometry.nu - dataset.geometry.mu
+        t0 = dataset.geometry.nu-dataset.geometry.mu
+        if not self.config.reduction.is_default('thetaRangeR'):
             # adjust range based on detector center
-            thetaRange = [ti+t0 for ti in self.reduction_config.thetaRangeR]
+            thetaRange = [ti+t0 for ti in self.config.reduction.thetaRangeR]
             proj.apply_theta_mask(thetaRange)
-        elif not self.reduction_config.is_default('thetaRange'):
-            proj.apply_theta_mask(self.reduction_config.thetaRange)
+        elif not self.config.reduction.is_default('thetaRange'):
+            proj.apply_theta_mask(self.config.reduction.thetaRange)
         else:
             thetaRange = [dataset.geometry.nu - dataset.geometry.mu - dataset.geometry.div/2,
                                               dataset.geometry.nu - dataset.geometry.mu + dataset.geometry.div/2]
             proj.apply_theta_mask(thetaRange)
+        for thi in self.config.reduction.thetaFilters:
+            # apply theta filters relative to angle on detector (issues with parts of the incoming divergence)
+            proj.apply_theta_filter((thi[0]+t0, thi[1]+t0))
 
-        proj.apply_lamda_mask(self.experiment_config.lambdaRange)
+        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)
 
-        if self.reduction_config.normalisationMethod == NormalisationMethod.over_illuminated:
+        if self.config.reduction.normalisationMethod == NormalisationMethod.over_illuminated:
             logging.debug('      assuming an overilluminated sample and correcting for the angle of incidence')
             proj.normalize_over_illuminated(self.norm)
-        elif self.reduction_config.normalisationMethod==NormalisationMethod.under_illuminated:
+        elif self.config.reduction.normalisationMethod==NormalisationMethod.under_illuminated:
             logging.debug('      assuming an underilluminated sample and ignoring the angle of incidence')
             proj.normalize_no_footprint(self.norm)
-        elif self.reduction_config.normalisationMethod==NormalisationMethod.direct_beam:
+        elif self.config.reduction.normalisationMethod==NormalisationMethod.direct_beam:
             logging.debug('      assuming direct beam for normalisation and ignoring the angle of incidence')
             proj.normalize_no_footprint(self.norm)
         else:

nicos_config.md

@@ -0,0 +1,43 @@
+EOS-Service
+===========
+
+EOS can be used as histogram service to send images to the Nicos instrument control software. 
+For that you need to run it on the amor instrument computer:
+
+```bash
+amor-nicos {-vv}
+```
+
+The instrument config in Nicos needs to configure a Kafka JustBinItImage instance 
+for each histogram that should be used:
+
+```python
+hist_yz = device('nicos_sinq.devices.just_bin_it.JustBinItImage',
+        description = 'Detector pixel histogram over all times',
+        hist_topic = 'AMOR_histograms_YZ',
+        data_topic = 'AMOR_detector',
+        command_topic = 'AMOR_histCommands',
+        brokers = ['linkafka01.psi.ch:9092'],
+        unit = 'evts',
+        hist_type = '2-D SANSLLB',
+        det_width = 446,
+        det_height = 64,
+        ),
+hist_tofz = device('nicos_sinq.devices.just_bin_it.JustBinItImage',
+        description = 'Detector time of flight vs. z-pixel histogram over all y-values',
+        hist_topic = 'AMOR_histograms_TofZ',
+        data_topic = 'AMOR_detector',
+        command_topic = 'AMOR_histCommands',
+        brokers = ['linkafka01.psi.ch:9092'],
+        unit = 'evts',
+        hist_type = '2-D SANSLLB',
+        det_width = 118,
+        det_height = 446,
+        ),
+```
+
+These images have then to be set in the detector configuration as _images_ items:
+
+```                      
+images=['hist_yz', 'hist_tofz'],
+```

requirements.txt

@@ -2,5 +2,7 @@ numpy
 h5py
 orsopy
 numba
+matplotlib
+tabulate
 backports.strenum; python_version<"3.11"
 backports.zoneinfo; python_version<"3.9"

setup.cfg

@@ -34,3 +34,6 @@ Homepage = "https://github.com/jochenstahn/amor"
 [options.entry_points]
 console_scripts =
     eos = eos.__main__:main
+    eosls = eos.ls:main
+    events2histogram = eos.e2h:main
+    amor-nicos = eos.nicos:main

tests/test_event_handling.py

@@ -0,0 +1,562 @@
+import os
+import numpy as np
+import logging
+
+from unittest import TestCase
+from datetime import datetime
+from copy import deepcopy
+
+from orsopy.fileio import Person, Experiment, Sample, InstrumentSettings, Value, ValueRange, Polarization
+
+from eos import const
+from eos.header import Header
+from eos.event_data_types import EVENT_BITMASKS, AmorGeometry, AmorTiming, AmorEventStream, \
+    EventDataAction, EventDatasetProtocol, PACKET_TYPE, PC_TYPE, PULSE_TYPE, EVENT_TYPE, append_fields
+from eos.event_handling import ApplyPhaseOffset, ApplyParameterOverwrites, CorrectChopperPhase, CorrectSeriesTime, \
+    AssociatePulseWithMonitor, FilterMonitorThreshold, FilterStrangeTimes, TofTimeCorrection, ApplyMask
+from eos.event_analysis import ExtractWalltime, MergeFrames, AnalyzePixelIDs, CalculateWavelength, CalculateQ, \
+    FilterQzRange, FilterByLog
+from eos.options import MonitorType, IncidentAngle, ExperimentConfig
+
+
+class MockEventData:
+    """
+    Simulated dataset to be used with event handling unit tests
+    """
+    geometry: AmorGeometry
+    timing: AmorTiming
+    data: AmorEventStream
+
+    def __init__(self):
+        self.geometry = AmorGeometry(mu=2.0, nu=1.0, kap=0.5, kad=0.0, div=1.5,
+                                     chopperSeparation=1000.0, detectorDistance=4000., chopperDetectorDistance=18842.)
+        self.timing = AmorTiming(
+                ch1TriggerPhase=-9.1, ch2TriggerPhase=6.75,
+                chopperPhase=0.17, chopperSpeed=500., tau=0.06
+                )
+        self.create_data()
+
+    def create_data(self):
+        # list of events, here with random time of fligh and pixel location
+        events = np.recarray((10000, ), dtype=EVENT_TYPE)
+        events.tof = np.random.uniform(low=0., high=0.12, size=events.shape)
+        events.pixelID = np.random.randint(0, 28671, size=events.shape)
+        events.mask = 0
+
+        # list of data packates containing previous events
+        packets = np.recarray((1000,), dtype=PACKET_TYPE)
+        packets.start_index = np.linspace(0, events.shape[0]-1, packets.shape[0], dtype=np.uint32)
+        packets.time = np.linspace(1700000000000000000, 1700000000000000000+3_600_000_000,
+                                   packets.shape[0], dtype=np.int64)
+
+        # chopper pulses within the measurement time
+        pulses = np.recarray((packets.shape[0],), dtype=PULSE_TYPE)
+        pulses.monitor = 1.0
+        pulses.time = packets.time
+
+        # proton current information with independent timing
+        proton_current = np.recarray((50,), dtype=PC_TYPE)
+        proton_current.current = 1500.0
+        proton_current[np.random.randint(0, proton_current.shape[0]-1, 10)] = 0. # random time with no current
+        proton_current.time = np.linspace(1700000000000000300, 1700000000000000000+3_600_000_000,
+                                   proton_current.shape[0], dtype=np.int64)
+
+        self.data = AmorEventStream(events, packets, pulses, proton_current)
+        self.orig_data = deepcopy(self.data)
+
+
+    def append(self, other):
+        raise NotImplementedError("Just for testing, no append")
+
+    def update_header(self, header:Header):
+        # update a header with the information read from file
+        header.owner = Person(name="test user", affiliation='PSI')
+        header.experiment = Experiment(title='test experiment', instrument='amor',
+                                       start_date=datetime.now(), probe="neutron")
+        header.sample = Sample(name='test sample')
+        header.measurement_instrument_settings = InstrumentSettings(incident_angle=Value(1.5, 'deg'),
+            wavelength = ValueRange(3.0, 12.5, 'angstrom'),
+            polarization = Polarization.unpolarized)
+
+    def update_info_from_logs(self):
+        RELEVANT_ITEMS = ['sample_temperature', 'sample_magnetic_field', 'polarization_config_label']
+        for key, log in self.data.device_logs.items():
+            if key not in RELEVANT_ITEMS:
+                continue
+            if log.value.dtype in [np.int8, np.int16, np.int32, np.int64]:
+                # for integer items (flags) report the most common one
+                value = np.bincount(log.value).argmax()
+                if logging.getLogger().getEffectiveLevel() <= logging.DEBUG \
+                        and np.unique(log.value).shape[0]>1:
+                    logging.debug(f'        filtered values for {key} not unique, '
+                                  f'has {np.unique(log.value).shape[0]} values')
+            else:
+                value = log.value.mean()
+            if key == 'polarization_config_label':
+                self.instrument_settings.polarization = Polarization(const.polarizationConfigs[value])
+            elif key == 'sample_temperature':
+                self.sample.sample_parameters['temperature'].magnitue = value
+            elif key == 'sample_magnetic_field':
+                self.sample.sample_parameters['magnetic_field'].magnitue = value
+
+
+class TestActionClass(TestCase):
+    @classmethod
+    def setUpClass(cls):
+        """
+        Create test classes to be used
+        """
+        class T1(EventDataAction):
+            def perform_action(self, event: EventDatasetProtocol):
+                event.data.events.mask += 1
+        class T2(EventDataAction):
+            def perform_action(self, event: EventDatasetProtocol):
+                event.data.events.mask += 2
+        class T4(EventDataAction):
+            def perform_action(self, event: EventDatasetProtocol):
+                event.data.events.mask += 4
+        cls.T1=T1; cls.T2=T2; cls.T4=T4
+
+        class H1(EventDataAction):
+            def perform_action(self, event: EventDatasetProtocol):
+                ...
+            def update_header(self, header:Header) ->None:
+                header.sample.name = 'h1'
+        class H2(EventDataAction):
+            def perform_action(self, event: EventDatasetProtocol):
+                ...
+            def update_header(self, header: Header) -> None:
+                header.sample.name = 'h2'
+        class HN(EventDataAction):
+            def __init__(self, n):
+                self._n = n
+            def perform_action(self, event: EventDatasetProtocol):
+                ...
+            def update_header(self, header: Header) -> None:
+                header.sample.name = self._n
+        cls.H1=H1; cls.H2=H2; cls.HN = HN
+
+    def setUp(self):
+        self.d = MockEventData()
+        self.header = Header()
+        self.d.update_header(self.header)
+
+    def test_individual(self):
+        t1 = self.T1()
+        t2 = self.T2()
+        t4 = self.T4()
+
+        np.testing.assert_array_equal(self.d.data.events.mask, 0)
+        t1.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, 1)
+        t2.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, 3)
+        t4.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, 7)
+
+    def test_header(self):
+        h1 = self.H1()
+        h2 = self.H2()
+        h3 = self.HN('h3')
+        h4 = self.HN('h4')
+
+        h1.update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h1')
+        h2.update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h2')
+        h3.update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h3')
+        h4.update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h4')
+
+    def test_combination(self):
+        t1 = self.T1()
+        t2 = self.T2()
+        t4 = self.T4()
+        t12 = t1 | t2
+        t24 = t2 | t4
+        t1224 = t1 | t2 | t2 | t4
+        t1224b = t12 | t24
+
+        np.testing.assert_array_equal(self.d.data.events.mask, 0)
+        t12.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, 3)
+        t24.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, 9)
+
+        t1224.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, 18)
+        t1224b.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, 27)
+
+
+    def test_combine_header(self):
+        h1 = self.H1()
+        h2 = self.H2()
+        h3 = self.HN('h3')
+        h4 = self.HN('h4')
+
+        (h1|h2).update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h2')
+        (h2|h1).update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h1')
+        (h3|h4).update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h4')
+        (h4|h3).update_header(self.header)
+        self.assertEqual(self.header.sample.name, 'h3')
+
+    def test_abstract_misssing(self):
+        with self.assertRaises(TypeError):
+            class E(EventDataAction):
+                ...
+            _ = E()
+
+    def test_hash(self):
+        """
+        Check that hashes of different actions are different but
+        instances of same action have same hash
+        """
+        t1 = self.T1()
+        t1b = self.T1()
+        t2 = self.T2()
+        t4 = self.T4()
+        h3 = self.HN('h3')
+        h3b = self.HN('h3')
+        h4 = self.HN('h4')
+
+        self.assertNotEqual(t1.action_hash(), t2.action_hash())
+        self.assertNotEqual(t2.action_hash(), t4.action_hash())
+        self.assertNotEqual(t1.action_hash(), t4.action_hash())
+        self.assertNotEqual(h3.action_hash(), h4.action_hash())
+        self.assertEqual(t1.action_hash(), t1b.action_hash())
+        self.assertEqual(h3.action_hash(), h3b.action_hash())
+
+
+class TestSimpleActions(TestCase):
+    def setUp(self):
+        self.d = MockEventData()
+        self.header = Header()
+        self.d.update_header(self.header)
+
+    def test_chopper_phase(self):
+        cp = CorrectChopperPhase()
+        cp.perform_action(self.d)
+        np.testing.assert_array_equal(
+                self.d.data.events.tof,
+                self.d.orig_data.events.tof+
+                self.d.timing.tau*(self.d.timing.ch1TriggerPhase-self.d.timing.chopperPhase/2)/180
+                )
+
+    def _extract_walltime(self):
+        # Extract wall time for events and orig copy
+        wt = ExtractWalltime()
+        d = self.d.data
+        self.d.data = self.d.orig_data
+        wt.perform_action(self.d)
+        self.d.data = d
+        wt.perform_action(self.d)
+
+    def test_extract_walltime(self):
+        self._extract_walltime()
+        # wallTime should be always a time present in the packet times
+        np.testing.assert_array_equal(np.isin(self.d.data.events.wallTime, self.d.data.packets.time), True)
+        # make sure extraction works on both original and copy
+        np.testing.assert_array_equal(self.d.data.events.wallTime, self.d.orig_data.events.wallTime)
+
+    def test_series_time(self):
+        corr = 100
+        ct = CorrectSeriesTime(corr)
+
+        with self.assertRaises(ValueError):
+            ct.perform_action(self.d)
+
+        self._extract_walltime()
+
+
+        ct.perform_action(self.d)
+        np.testing.assert_array_equal(
+                self.d.data.pulses.time,
+                self.d.orig_data.pulses.time-corr
+                )
+        np.testing.assert_array_equal(
+                self.d.data.events.wallTime,
+                self.d.orig_data.events.wallTime-corr
+                )
+        np.testing.assert_array_equal(
+                self.d.data.proton_current.time,
+                self.d.orig_data.proton_current.time-corr
+                )
+
+    def test_associate_monitor(self):
+        amPC = AssociatePulseWithMonitor(MonitorType.proton_charge)
+        amT = AssociatePulseWithMonitor(MonitorType.time)
+        amN = AssociatePulseWithMonitor(MonitorType.neutron_monitor)
+
+        self.d.data.pulses.monitor = 13
+        amN.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.pulses.monitor, 1)
+
+        self.d.data.pulses.monitor = 13
+        amT.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.pulses.monitor, np.float32(2*self.d.timing.tau))
+
+        self.d.data.pulses.monitor = 13
+        amPC.perform_action(self.d)
+        pcm = self.d.data.proton_current.current *2*self.d.timing.tau*1e-3
+        np.testing.assert_array_equal(np.isin(self.d.data.pulses.monitor, pcm), True)
+
+    def test_filter_monitor_threashold(self):
+        amPC = AssociatePulseWithMonitor(MonitorType.proton_charge)
+        fmt = amPC | FilterMonitorThreshold(1000.)
+        fma = amPC | FilterMonitorThreshold(2000.)
+        fm0 = amPC | FilterMonitorThreshold(-1.0)
+
+        with self.assertRaises(ValueError):
+            fmt.perform_action(self.d)
+
+        self._extract_walltime()
+        fm0.perform_action(self.d)
+        self.assertEqual(self.d.data.events.mask.sum(), 0)
+        fmt.perform_action(self.d)
+        # calculate, which events should have 0 monitor
+        zero_times = self.d.data.pulses.time[self.d.data.pulses.monitor==0]
+        zero_sum = np.isin(self.d.data.events.wallTime, zero_times).sum()
+        self.assertEqual(self.d.data.events.mask.sum(), zero_sum*EVENT_BITMASKS['MonitorThreshold'])
+        # filter all events
+        self.d.data.events.mask = 0
+        fma.perform_action(self.d)
+        self.assertEqual(self.d.data.events.mask.sum(), self.d.data.events.shape[0]*EVENT_BITMASKS['MonitorThreshold'])
+
+    def test_filter_strage_times(self):
+        st = FilterStrangeTimes()
+
+        st.perform_action(self.d)
+        self.assertEqual(self.d.data.events.mask.sum(), 0)
+
+        # half events should be strange times (outside of ToF frame)
+        self.d.data.events.tof += self.d.timing.tau
+        st.perform_action(self.d)
+        self.assertEqual(self.d.data.events.mask.sum(),
+                         (self.d.data.events.tof>2*self.d.timing.tau).sum()*EVENT_BITMASKS['StrangeTimes'])
+
+    def test_apply_phase_offset(self):
+        action = ApplyPhaseOffset(12.5)
+        action.perform_action(self.d)
+        self.assertEqual(self.d.timing.ch1TriggerPhase, 12.5)
+
+    def test_apply_parameter_overwrites(self):
+        action = ApplyParameterOverwrites(ExperimentConfig(muOffset=0.25, mu=3.5, nu=4.5))
+        action.perform_action(self.d)
+        self.assertEqual(self.d.geometry.mu, 3.5)
+        self.assertEqual(self.d.geometry.nu, 4.5)
+
+        action = ApplyParameterOverwrites(ExperimentConfig(muOffset=0.25))
+        action.perform_action(self.d)
+        self.assertEqual(self.d.geometry.mu, 3.75)
+
+        action = ApplyParameterOverwrites(ExperimentConfig(sampleModel='air | Si | Fe'))
+        action.update_header(self.header)
+        self.assertIsNotNone(self.header.sample.model)
+
+    def test_apply_sample_model_file(self):
+        if os.path.isfile('test.yaml'):
+            os.remove('test.yaml')
+        action = ApplyParameterOverwrites(ExperimentConfig(sampleModel='test.yaml'))
+        action.update_header(self.header)
+        self.assertIsNone(self.header.sample.model)
+
+        with open('test.yaml', 'w') as fh:
+            fh.write("""stack: air | Si | Fe""")
+
+        try:
+            action = ApplyParameterOverwrites(ExperimentConfig(sampleModel='test.yaml'))
+            action.update_header(self.header)
+            self.assertEqual(self.header.sample.model.stack, 'air | Si | Fe')
+        finally:
+            os.remove('test.yaml')
+
+    def test_tof_time_correction(self):
+        action = TofTimeCorrection()
+        with self.assertRaises(ValueError):
+            action.perform_action(self.d)
+
+        new_events = append_fields(self.d.data.events, [('delta', np.float64)])
+        new_events.delta = 10.0
+        self.d.data.events = new_events
+        tof_before = self.d.data.events.tof.copy()
+        action.perform_action(self.d)
+        np.testing.assert_allclose(
+            self.d.data.events.tof,
+            tof_before - (10.0 / 180.0) * self.d.timing.tau
+        )
+
+        self.d.create_data()
+        new_events = append_fields(self.d.data.events, [('delta', np.float64)])
+        new_events.delta = 10.0
+        self.d.data.events = new_events
+        tof_before = self.d.data.events.tof.copy()
+        action = TofTimeCorrection(correct_chopper_opening=False)
+        action.perform_action(self.d)
+        np.testing.assert_allclose(
+            self.d.data.events.tof,
+            tof_before - (self.d.geometry.kad / 180.0) * self.d.timing.tau
+        )
+
+    def test_apply_mask(self):
+        self.d.data.events = self.d.data.events[:6].copy()
+        self.d.data.events.mask[:] = [0, 1, 2, 3, 4, 5]
+
+        action = ApplyMask()
+        action.perform_action(self.d)
+        self.assertEqual(self.d.data.events.shape[0], 1)
+        self.assertEqual(self.d.data.events.mask[0], 0)
+
+        self.d.create_data()
+        self.d.data.events = self.d.data.events[:6].copy()
+        self.d.data.events.mask[:] = [0, 1, 2, 3, 4, 5]
+        action = ApplyMask(bitmask_filter=EVENT_BITMASKS['MonitorThreshold'])
+        action.perform_action(self.d)
+        np.testing.assert_array_equal(self.d.data.events.mask, np.array([0, EVENT_BITMASKS['MonitorThreshold']],
+                                                                        dtype=np.int32))
+
+    def test_merge_frames(self):
+        action = MergeFrames(lamdaCut=0.0)
+        action.perform_action(self.d)
+        self.assertEqual(self.d.data.events.tof.shape, self.d.orig_data.events.tof.shape)
+        np.testing.assert_array_compare(lambda x,y: x<=y, self.d.data.events.tof, self.d.orig_data.events.tof)
+        self.assertTrue((-self.d.timing.tau<=self.d.data.events.tof).all())
+        np.testing.assert_array_less(self.d.data.events.tof, self.d.timing.tau)
+
+        action = MergeFrames(lamdaCut=2.0)
+        self.d.data.events.tof = self.d.orig_data.events.tof[:]
+        action.perform_action(self.d)
+        tofCut = 2.0*self.d.geometry.chopperDetectorDistance/const.hdm*1e-13
+        self.assertTrue((tofCut-self.d.timing.tau<=self.d.data.events.tof).all())
+        self.assertTrue((self.d.data.events.tof<=tofCut+self.d.timing.tau).all())
+
+    def test_analyze_pixel_ids(self):
+        action = AnalyzePixelIDs((1000, 1001))
+        action.perform_action(self.d)
+        self.assertIn('detZ', self.d.data.events.dtype.names)
+        self.assertIn('detXdist', self.d.data.events.dtype.names)
+        self.assertIn('delta', self.d.data.events.dtype.names)
+        self.assertEqual(
+            np.bitwise_and(self.d.data.events.mask, EVENT_BITMASKS['yRange']).astype(bool).sum(),
+            self.d.data.events.shape[0]
+        )
+        # TODO: maybe add a test actually checking correct detector-id resolution
+
+    def test_calculate_wavelength(self):
+        action = CalculateWavelength((3.0, 5.0))
+        with self.assertRaises(ValueError):
+            action.perform_action(self.d)
+
+        new_events = append_fields(self.d.data.events, [('detXdist', np.float64)])
+        new_events.detXdist = 0.0
+        self.d.data.events = new_events
+        action.perform_action(self.d)
+        self.assertIn('lamda', self.d.data.events.dtype.names)
+        flt = self.d.data.events.mask!=EVENT_BITMASKS['LamdaRange']
+        # check all wavelength in range not filtered
+        np.testing.assert_array_less(self.d.data.events.lamda[flt], 5.0)
+        np.testing.assert_array_less(3.0, self.d.data.events.lamda[flt])
+        # check all wavelength out of range filtered
+        flt = self.d.data.events.mask==EVENT_BITMASKS['LamdaRange']
+        self.assertTrue(((self.d.data.events.lamda[flt]<3.0)|(self.d.data.events.lamda[flt]>5.0)).all())
+
+    def test_calculate_q(self):
+        action = CalculateQ(IncidentAngle.alphaF)
+        with self.assertRaises(ValueError):
+            action.perform_action(self.d)
+
+        # TODO: add checks for actual resulting values
+
+        new_events = append_fields(self.d.data.events, [('lamda', np.float64), ('delta', np.float64)])
+        new_events.lamda = 5.0
+        new_events.delta = 0.0
+        self.d.data.events = new_events
+        action.perform_action(self.d)
+        self.assertIn('qz', self.d.data.events.dtype.names)
+        self.assertNotIn('qx', self.d.data.events.dtype.names)
+        action.update_header(self.header)
+        self.assertEqual(self.header.measurement_scheme, 'angle- and energy-dispersive')
+
+        self.d.create_data()
+        new_events = append_fields(self.d.data.events, [('lamda', np.float64), ('delta', np.float64)])
+        new_events.lamda = 5.0
+        new_events.delta = 0.0
+        self.d.data.events = new_events
+        action = CalculateQ(IncidentAngle.mu)
+        action.perform_action(self.d)
+        self.assertIn('qz', self.d.data.events.dtype.names)
+        self.assertIn('qx', self.d.data.events.dtype.names)
+        action.update_header(self.header)
+        self.assertEqual(self.header.measurement_scheme, 'energy-dispersive')
+
+        self.d.create_data()
+        new_events = append_fields(self.d.data.events, [('lamda', np.float64), ('delta', np.float64)])
+        new_events.lamda = 5.0
+        new_events.delta = 0.0
+        self.d.data.events = new_events
+        action = CalculateQ(IncidentAngle.nu)
+        action.perform_action(self.d)
+        self.assertIn('qz', self.d.data.events.dtype.names)
+        self.assertNotIn('qx', self.d.data.events.dtype.names)
+        action.update_header(self.header)
+        self.assertEqual(self.header.measurement_scheme, 'energy-dispersive')
+
+    def test_filter_qz_range(self):
+        action = FilterQzRange((0.1, 0.2))
+        with self.assertRaises(ValueError):
+            action.perform_action(self.d)
+
+        self.d.data.events = self.d.data.events[:5].copy()
+        new_events = append_fields(self.d.data.events, [('qz', np.float64)])
+        new_events.qz = np.array([0.05, 0.1, 0.15, 0.2, 0.25])
+        self.d.data.events = new_events
+        action.perform_action(self.d)
+        np.testing.assert_array_equal(
+            self.d.data.events.mask,
+            np.array([1, 0, 0, 0, 1], dtype=np.int32) * EVENT_BITMASKS['qRange']
+        )
+
+    def test_filter_by_log(self):
+        action = FilterByLog("test_log==0") | ApplyMask()
+        class LogWarnError(Exception):
+            ...
+        def warn_raise(*args, **kwargs):
+            raise LogWarnError()
+        _orig_warn = logging.warning
+        try:
+            logging.warning = warn_raise
+            with self.assertRaises(LogWarnError):
+                action.perform_action(self.d)
+        finally:
+            logging.warning = _orig_warn
+
+        self._extract_walltime()
+
+        test_log = np.recarray(shape=(2,), dtype=np.dtype([('value', np.int32),
+                                                           ('time', np.int64)]))
+        test_log.time = [-5, self.d.data.pulses.time[100]+123]
+        test_log.value = [0, 1]
+        self.d.data.device_logs['test_log'] = test_log
+        action.perform_action(self.d)
+        self.assertEqual(self.d.data.pulses.shape[0], 101)
+
+    def test_filter_by_log_switchpulse(self):
+        action = FilterByLog("!test_log==0") | ApplyMask()
+        self._extract_walltime()
+
+        test_log = np.recarray(shape=(2,), dtype=np.dtype([('value', np.int32),
+                                                           ('time', np.int64)]))
+        test_log.time = [-5, self.d.data.pulses.time[100]+123]
+        test_log.value = [0, 1]
+        self.d.data.device_logs['test_log'] = test_log
+        self.d.data.device_logs['check_log'] = test_log.copy()
+        action.perform_action(self.d)
+        self.assertEqual(self.d.data.pulses.shape[0], 100)
+        np.testing.assert_array_equal(
+                self.d.data.device_logs['test_log'],
+                self.d.data.device_logs['check_log'],
+                )

tests/test_full_analysis.py

@@ -1,8 +1,10 @@
 import os
 import cProfile
+import numpy as np
 from unittest import TestCase
 from dataclasses import fields, MISSING
-from eos import options, reduction, logconfig
+from eos import options, reduction_reflectivity, logconfig
+from orsopy import fileio
 
 logconfig.setup_logging()
 logconfig.update_loglevel(1)
@@ -14,7 +16,7 @@ class FullAmorTest(TestCase):
     def setUpClass(cls):
         # generate map for option defaults
         cls._field_defaults = {}
-        for opt in [options.ExperimentConfig, options.ReductionConfig, options.OutputConfig]:
+        for opt in [options.ExperimentConfig, options.ReflectivityReductionConfig, options.ReflectivityOutputConfig]:
             defaults = {}
             for field in fields(opt):
                 if field.default not in [None, MISSING]:
@@ -38,82 +40,125 @@ class FullAmorTest(TestCase):
     def tearDown(self):
         self.pr.disable()
         for fi in ['test_results/test.Rqz.ort', 'test_results/5952.norm']:
-            try:
-                os.unlink(fi)
-            except FileNotFoundError:
-                pass
+           try:
+               os.unlink(fi)
+           except FileNotFoundError:
+               pass
 
     def test_time_slicing(self):
         experiment_config = options.ExperimentConfig(
-                chopperSpeed=self._field_defaults['ExperimentConfig']['chopperSpeed'],
                 chopperPhase=-13.5,
                 chopperPhaseOffset=-5,
-                monitorType=self._field_defaults['ExperimentConfig']['monitorType'],
-                lowCurrentThreshold=self._field_defaults['ExperimentConfig']['lowCurrentThreshold'],
                 yRange=(18, 48),
                 lambdaRange=(3., 11.5),
-                incidentAngle=self._field_defaults['ExperimentConfig']['incidentAngle'],
                 mu=0,
                 nu=0,
                 muOffset=0.0,
                 sampleModel='air | 10 H2O | D2O'
                 )
-        reduction_config = options.ReductionConfig(
-                normalisationMethod=self._field_defaults['ReductionConfig']['normalisationMethod'],
+        reduction_config = options.ReflectivityReductionConfig(
                 qResolution=0.01,
-                qzRange=self._field_defaults['ReductionConfig']['qzRange'],
+                qzRange=(0.01, 0.15),
                 thetaRange=(-0.75, 0.75),
                 fileIdentifier=["6003-6005"],
                 scale=[1],
                 normalisationFileIdentifier=[],
                 timeSlize=[300.0]
                 )
-        output_config = options.OutputConfig(
+        output_config = options.ReflectivityOutputConfig(
                 outputFormats=[options.OutputFomatOption.Rqz_ort],
                 outputName='test',
                 outputPath='test_results',
                 )
-        config=options.EOSConfig(self.reader_config, experiment_config, reduction_config, output_config)
+        config=options.ReflectivityConfig(self.reader_config, experiment_config, reduction_config, output_config)
         # run three times to get similar timing to noslicing runs
-        reducer = reduction.AmorReduction(config)
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
         reducer.reduce()
-        reducer = reduction.AmorReduction(config)
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
         reducer.reduce()
-        reducer = reduction.AmorReduction(config)
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
         reducer.reduce()
 
     def test_noslicing(self):
         experiment_config = options.ExperimentConfig(
-                chopperSpeed=self._field_defaults['ExperimentConfig']['chopperSpeed'],
                 chopperPhase=-13.5,
                 chopperPhaseOffset=-5,
-                monitorType=self._field_defaults['ExperimentConfig']['monitorType'],
-                lowCurrentThreshold=self._field_defaults['ExperimentConfig']['lowCurrentThreshold'],
                 yRange=(18, 48),
                 lambdaRange=(3., 11.5),
-                incidentAngle=self._field_defaults['ExperimentConfig']['incidentAngle'],
                 mu=0,
                 nu=0,
                 muOffset=0.0,
                 )
-        reduction_config = options.ReductionConfig(
-                normalisationMethod=self._field_defaults['ReductionConfig']['normalisationMethod'],
+        reduction_config = options.ReflectivityReductionConfig(
                 qResolution=0.01,
-                qzRange=self._field_defaults['ReductionConfig']['qzRange'],
                 thetaRange=(-0.75, 0.75),
                 fileIdentifier=["6003", "6004", "6005"],
                 scale=[1],
                 normalisationFileIdentifier=["5952"],
                 autoscale=(0.0, 0.05),
                 )
-        output_config = options.OutputConfig(
+        output_config = options.ReflectivityOutputConfig(
                 outputFormats=[options.OutputFomatOption.Rqz_ort],
                 outputName='test',
                 outputPath='test_results',
                 )
-        config=options.EOSConfig(self.reader_config, experiment_config, reduction_config, output_config)
-        reducer = reduction.AmorReduction(config)
+        config=options.ReflectivityConfig(self.reader_config, experiment_config, reduction_config, output_config)
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
         reducer.reduce()
         # run second time to reuse norm file
-        reducer = reduction.AmorReduction(config)
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
         reducer.reduce()
+
+    def test_eventfilter(self):
+        self.reader_config.year = 2026
+        experiment_config = options.ExperimentConfig()
+        reduction_config = options.ReflectivityReductionConfig(fileIdentifier=["826"],
+                                                               logfilter=['polarization_config_label==2'])
+        output_config = options.ReflectivityOutputConfig(
+                outputFormats=[options.OutputFomatOption.Rqz_ort],
+                outputName='test',
+                outputPath='test_results',
+                )
+        config=options.ReflectivityConfig(self.reader_config, experiment_config, reduction_config, output_config)
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
+        reducer.reduce()
+        espin_up = reducer.dataset.data.events.shape[0]
+
+        reduction_config.logfilter = ['polarization_config_label==3']
+        output_config.append = True
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
+        reducer.reduce()
+        espin_down = reducer.dataset.data.events.shape[0]
+        # measurement should have about 2x as many counts in spin_down
+        self.assertAlmostEqual(espin_down/espin_up, 2., 2)
+
+        # perform the same filter but remove pulses during which the switch occured
+        reduction_config.logfilter = ['!polarization_config_label==3']
+        output_config.append = True
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
+        reducer.reduce()
+        espin_down2 = reducer.dataset.data.events.shape[0]
+        # measurement should have about 2x as many counts in spin_down
+        self.assertLess(espin_down2, espin_down)
+
+    def test_polsplitting(self):
+        self.reader_config.year = 2026
+        experiment_config = options.ExperimentConfig()
+        reduction_config = options.ReflectivityReductionConfig(fileIdentifier=["826"])
+        output_config = options.ReflectivityOutputConfig(
+                outputFormats=[options.OutputFomatOption.Rqz_ort],
+                outputName='test',
+                outputPath='test_results',
+                )
+        config=options.ReflectivityConfig(self.reader_config, experiment_config, reduction_config, output_config)
+        reducer = reduction_reflectivity.ReflectivityReduction(config)
+        reducer.reduce()
+
+        results = fileio.load_orso(os.path.join(output_config.outputPath, output_config.outputName+'.Rqz.ort'))
+        self.assertEqual(len(results), 2)
+        self.assertEqual(results[0].info.data_source.measurement.instrument_settings.polarization, 'po')
+        self.assertEqual(results[1].info.data_source.measurement.instrument_settings.polarization, 'mo')
+        espin_up = np.nansum(results[0].data[:,1])
+        espin_down = np.nansum(results[1].data[:,1])
+        # the total intensity should be around equal as events are doubled and monitor counts are doubled
+        self.assertAlmostEqual(espin_down/espin_up, 1., 2)

update.md

@@ -0,0 +1,16 @@
+Make new release
+================
+
+- Update revision in `eos/__init__.py`
+- Commit changes `git commit -a -m "your message here"`
+- Tag version `git tag v3.x.y`
+- Push changes `git push` and `git push --tags`
+- This should trigger the **Release** action on GitHub that builds a new version and uploads it to PyPI.
+
+
+Update on AMOR
+==============
+
+- Login via SSH using the **amor** user.
+- Activate eos virtual environment `source /home/software/virtualenv/eosenv/bin/activate`
+- Update eos packge `pip install --upgrade amor-eos`