Compare commits
65 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
| f55d840b8e | |||
| 327d6a0ff8 | |||
| 22943bc513 | |||
|
|
01873d60db | ||
| 654251a194 | |||
| 9e30970d9b | |||
| 7b6f2045cc | |||
| e3d15a1142 | |||
|
|
98244327eb | ||
| 61ba1cf084 | |||
| 54f2169888 | |||
| 24cc0a4287 | |||
| a3d342e5d7 | |||
| c92e4e4d17 | |||
| 4d74237669 | |||
| 6144200551 | |||
| f2c681ffba | |||
| a56f29191d | |||
| 6597ca22d1 | |||
| ce31dcb190 | |||
| 41c5218413 | |||
| 8abd977656 | |||
| 432d85c9b3 | |||
| c222f42a89 | |||
| fbced41b9f | |||
| dfb5aa319f | |||
| c073fae269 | |||
| 447b81d09d | |||
| 003a8c04ce | |||
| 4cd25dec3d | |||
| b38eeb8eb3 | |||
| d722228079 | |||
| 493c2bf802 | |||
| c429429dd2 | |||
| 0550e725e8 | |||
| 5c8b9a8cd6 | |||
| 4e5d085ad7 | |||
| 2d8d1c66c6 | |||
| 14ea89e243 | |||
| 322c00ca54 | |||
| 0b7a56628f | |||
| 77641412ab | |||
| fc1bd66c3d | |||
| 95a1ffade4 | |||
| 4ee1cf7ea7 | |||
| c90bdd3316 | |||
| ac9babb9ad | |||
| aa5b540aed | |||
| 86d676ccc8 | |||
| f8daa93c48 | |||
| 4c2e344d78 | |||
| 469e1c0ab0 | |||
| a947dfd398 | |||
| 4e2269bdae | |||
| 3f93ec2017 | |||
| 9d788da2f1 | |||
| f16feac748 | |||
| 709a0c5392 | |||
| b71b72c6a3 | |||
| 2972ffd5d8 | |||
| 31201795b0 | |||
| 99af021b3e | |||
| aacbe3ed6f | |||
| 6c0c2fcab8 | |||
| db5713ab56 |
30
.github/workflows/release.yml
vendored
30
.github/workflows/release.yml
vendored
@@ -22,7 +22,36 @@ on:
|
||||
- all_incl_release
|
||||
|
||||
jobs:
|
||||
test:
|
||||
|
||||
runs-on: ubuntu-22.04
|
||||
strategy:
|
||||
matrix:
|
||||
python-version: ['3.8', '3.9', '3.10', '3.11', '3.12', '3.13']
|
||||
fail-fast: false
|
||||
|
||||
steps:
|
||||
- uses: actions/checkout@v4
|
||||
with:
|
||||
lfs: 'true'
|
||||
- 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:
|
||||
@@ -55,6 +84,7 @@ jobs:
|
||||
skip-existing: true
|
||||
|
||||
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)) }}
|
||||
|
||||
|
||||
76
.github/workflows/unit_tests.yml
vendored
76
.github/workflows/unit_tests.yml
vendored
@@ -1,38 +1,38 @@
|
||||
name: Unit Testing
|
||||
|
||||
on:
|
||||
push:
|
||||
branches:
|
||||
- main
|
||||
pull_request:
|
||||
|
||||
# Allows you to run this workflow manually from the Actions tab
|
||||
workflow_dispatch:
|
||||
|
||||
jobs:
|
||||
build:
|
||||
|
||||
runs-on: ubuntu-22.04
|
||||
strategy:
|
||||
matrix:
|
||||
python-version: ['3.8', '3.9', '3.10', '3.11', '3.12']
|
||||
fail-fast: false
|
||||
|
||||
steps:
|
||||
- uses: actions/checkout@v4
|
||||
with:
|
||||
lfs: 'true'
|
||||
- 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
|
||||
name: Unit Testing
|
||||
|
||||
on:
|
||||
push:
|
||||
branches:
|
||||
- main
|
||||
pull_request:
|
||||
|
||||
# Allows you to run this workflow manually from the Actions tab
|
||||
workflow_dispatch:
|
||||
|
||||
jobs:
|
||||
test:
|
||||
|
||||
runs-on: ubuntu-22.04
|
||||
strategy:
|
||||
matrix:
|
||||
python-version: ['3.8', '3.9', '3.10', '3.11', '3.12']
|
||||
fail-fast: false
|
||||
|
||||
steps:
|
||||
- uses: actions/checkout@v4
|
||||
with:
|
||||
lfs: 'true'
|
||||
- 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
|
||||
|
||||
@@ -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.1.1'
|
||||
__date__ = '2026-02-26'
|
||||
|
||||
@@ -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()
|
||||
|
||||
|
||||
@@ -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):
|
||||
"""
|
||||
Process command line argument.
|
||||
The type of the default values is used for conversion and validation.
|
||||
|
||||
42
eos/e2h.py
Normal file
42
eos/e2h.py
Normal file
@@ -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()
|
||||
@@ -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)
|
||||
|
||||
@@ -31,7 +31,7 @@ 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)])
|
||||
|
||||
|
||||
@@ -77,49 +77,28 @@ class CorrectSeriesTime(EventDataAction):
|
||||
|
||||
|
||||
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 +113,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 +149,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
|
||||
|
||||
@@ -29,94 +29,77 @@ except ImportError:
|
||||
AMOR_LOCAL_TIMEZONE = zoneinfo.ZoneInfo(key='Europe/Zurich')
|
||||
|
||||
if platform.node().startswith('amor'):
|
||||
NICOS_CACHE_DIR = '/home/amor/nicosdata/amor/cache/'
|
||||
GREP = '/usr/bin/grep "%s"'
|
||||
NICOS_CACHE_DIR = '/home/data/nicosdata/cache/'
|
||||
GREP = '/usr/bin/grep "value"'
|
||||
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
|
||||
|
||||
eventStartTime: np.int64
|
||||
|
||||
def __init__(self, fileName:Union[str, h5py.File, BinaryIO], first_index:int=0, max_events:int=100_000_000):
|
||||
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.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=''):
|
||||
try:
|
||||
return dtype(self.hdf[f'/entry1/Amor/{key}'][0])
|
||||
except(KeyError, IndexError):
|
||||
if NICOS_CACHE_DIR:
|
||||
try:
|
||||
logging.warning(f" using parameter {nicos_key} from nicos cache")
|
||||
logging.info(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]
|
||||
call = f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year_date}{suffix}'
|
||||
value = str(subprocess.getoutput(call)).split('\t')[-1]
|
||||
return dtype(value)
|
||||
except Exception:
|
||||
logging.error("Couldn't get value from nicos cache", exc_info=True)
|
||||
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 get_hdf_single_entry(self, path):
|
||||
if not np.shape(self.hdf['entry1/title']):
|
||||
return self.hdf[path][()].decode('utf-8')
|
||||
else:
|
||||
# format until 2025
|
||||
return self.hdf[path][0].decode('utf-8')
|
||||
|
||||
def read_header_info(self):
|
||||
# read general information and first data set
|
||||
title = self.hdf['entry1/title'][0].decode('utf-8')
|
||||
proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8')
|
||||
user_name = self.hdf['entry1/user/name'][0].decode('utf-8')
|
||||
title = self.get_hdf_single_entry('entry1/title')
|
||||
proposal_id = self.get_hdf_single_entry('entry1/proposal_id')
|
||||
user_name = self.get_hdf_single_entry('entry1/user/name')
|
||||
user_affiliation = 'unknown'
|
||||
user_email = self.hdf['entry1/user/email'][0].decode('utf-8')
|
||||
user_email = self.get_hdf_single_entry('entry1/user/email')
|
||||
user_orcid = None
|
||||
sampleName = self.hdf['entry1/sample/name'][0].decode('utf-8')
|
||||
model = self.hdf['entry1/sample/model'][0].decode('utf-8')
|
||||
if 'stack:' in model:
|
||||
sampleName = self.get_hdf_single_entry('entry1/sample/name')
|
||||
instrumentName = 'Amor'
|
||||
source = self.get_hdf_single_entry('entry1/Amor/source/name')
|
||||
sourceProbe = 'neutron'
|
||||
model = self.get_hdf_single_entry('entry1/sample/model')
|
||||
if 'stack' in model:
|
||||
import yaml
|
||||
model = yaml.safe_load(model)
|
||||
else:
|
||||
model = dict(stack=model)
|
||||
instrumentName = 'Amor'
|
||||
source = self.hdf['entry1/Amor/source/name'][0].decode('utf-8')
|
||||
sourceProbe = 'neutron'
|
||||
start_time = self.hdf['entry1/start_time'][0].decode('utf-8')
|
||||
start_time = self.get_hdf_single_entry('entry1/start_time')
|
||||
# extract start time as unix time, adding UTC offset of 1h to time string
|
||||
start_date = datetime.fromisoformat(start_time)
|
||||
self.fileDate = start_date.replace(tzinfo=AMOR_LOCAL_TIMEZONE)
|
||||
@@ -137,11 +120,23 @@ 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 \
|
||||
and float(self.hdf['entry1/sample/temperature/average_value'][0])>0:
|
||||
self.sample.sample_parameters['temperature'] = fileio.Value(
|
||||
float(self.hdf['entry1/sample/temperature/average_value'][0]), unit='K')
|
||||
if self.hdf['entry1/sample'].get('magnetic_field', None) is not None:
|
||||
self.sample.sample_parameters['magnetic_field'] = fileio.Value(
|
||||
float(self.hdf['entry1/sample/magnetic_field/average_value'][0]), unit='T')
|
||||
|
||||
def read_instrument_configuration(self):
|
||||
chopperSeparation = float(np.take(self.hdf['entry1/Amor/chopper/pair_separation'], 0))
|
||||
@@ -153,7 +148,7 @@ class AmorEventData:
|
||||
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)
|
||||
kad = self._replace_if_missing('instrument_control_parameters/kappa_offset', 'kappa_offset', 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)
|
||||
@@ -177,7 +172,7 @@ 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)
|
||||
polarizationConfigLabel = self._replace_if_missing('polarization/configuration/average_value', 'polarization_config_label', int, suffix='/*')
|
||||
polarizationConfig = fileio.Polarization(polarizationConfigs[polarizationConfigLabel])
|
||||
logging.debug(f' polarization configuration: {polarizationConfig} (index {polarizationConfigLabel})')
|
||||
|
||||
@@ -187,9 +182,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 +220,52 @@ 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
|
||||
self.read_event_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 read_event_stream(self):
|
||||
"""
|
||||
Read the actual event data from file. If file is too large, find event index from packets
|
||||
@@ -230,12 +273,13 @@ 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:]
|
||||
|
||||
nevts = self.hdf['/entry1/Amor/detector/data/event_time_offset'].shape[0]
|
||||
@@ -256,15 +300,15 @@ 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_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)
|
||||
@@ -281,17 +325,17 @@ class AmorEventData:
|
||||
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])]
|
||||
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.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):
|
||||
|
||||
@@ -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
|
||||
#-------------------------------------------------------------------------------------------------
|
||||
|
||||
@@ -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
|
||||
|
||||
165
eos/kafka_events.py
Normal file
165
eos/kafka_events.py
Normal file
@@ -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)
|
||||
283
eos/kafka_serializer.py
Normal file
283
eos/kafka_serializer.py
Normal file
@@ -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 == self._active_histogram_yz:
|
||||
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()
|
||||
46
eos/nicos.py
Normal file
46
eos/nicos.py
Normal file
@@ -0,0 +1,46 @@
|
||||
"""
|
||||
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],
|
||||
'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()
|
||||
@@ -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
|
||||
|
||||
340
eos/options.py
340
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,33 +471,6 @@ 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))
|
||||
|
||||
|
||||
class OutputFomatOption(StrEnum):
|
||||
Rqz_ort = "Rqz.ort"
|
||||
@@ -440,9 +489,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={
|
||||
@@ -511,11 +561,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 +578,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
|
||||
|
||||
@@ -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}')
|
||||
|
||||
|
||||
@@ -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,25 @@ 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)
|
||||
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 +249,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 +257,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 +289,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):
|
||||
# ...
|
||||
318
eos/reduction_e2h.py
Normal file
318
eos/reduction_e2h.py
Normal file
@@ -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=1000)
|
||||
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)
|
||||
128
eos/reduction_kafka.py
Normal file
128
eos/reduction_kafka.py
Normal file
@@ -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)
|
||||
@@ -8,27 +8,21 @@ from orsopy import fileio
|
||||
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 +31,66 @@ 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)
|
||||
self.dataevent_actions |= ea.FilterQzRange(self.config.reduction.qzRange)
|
||||
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 +100,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()
|
||||
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,11 +159,11 @@ 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 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)
|
||||
@@ -162,26 +171,26 @@ class AmorReduction:
|
||||
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]}')
|
||||
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'
|
||||
|
||||
# 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 +205,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 +252,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):
|
||||
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 +277,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:
|
||||
@@ -303,11 +312,11 @@ class AmorReduction:
|
||||
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 +324,19 @@ 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)
|
||||
|
||||
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 +349,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 +373,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:
|
||||
1056
events2histogram.py
1056
events2histogram.py
File diff suppressed because it is too large
Load Diff
43
nicos_config.md
Normal file
43
nicos_config.md
Normal file
@@ -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'],
|
||||
```
|
||||
@@ -2,5 +2,6 @@ numpy
|
||||
h5py
|
||||
orsopy
|
||||
numba
|
||||
matplotlib
|
||||
backports.strenum; python_version<"3.11"
|
||||
backports.zoneinfo; python_version<"3.9"
|
||||
|
||||
@@ -34,3 +34,5 @@ Homepage = "https://github.com/jochenstahn/amor"
|
||||
[options.entry_points]
|
||||
console_scripts =
|
||||
eos = eos.__main__:main
|
||||
events2histogram = eos.e2h:main
|
||||
amor-nicos = eos.nicos:main
|
||||
|
||||
@@ -2,7 +2,7 @@ import os
|
||||
import cProfile
|
||||
from unittest import TestCase
|
||||
from dataclasses import fields, MISSING
|
||||
from eos import options, reduction, logconfig
|
||||
from eos import options, reduction_reflectivity, logconfig
|
||||
|
||||
logconfig.setup_logging()
|
||||
logconfig.update_loglevel(1)
|
||||
@@ -14,7 +14,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]:
|
||||
@@ -58,28 +58,28 @@ class FullAmorTest(TestCase):
|
||||
muOffset=0.0,
|
||||
sampleModel='air | 10 H2O | D2O'
|
||||
)
|
||||
reduction_config = options.ReductionConfig(
|
||||
normalisationMethod=self._field_defaults['ReductionConfig']['normalisationMethod'],
|
||||
reduction_config = options.ReflectivityReductionConfig(
|
||||
normalisationMethod=self._field_defaults['ReflectivityReductionConfig']['normalisationMethod'],
|
||||
qResolution=0.01,
|
||||
qzRange=self._field_defaults['ReductionConfig']['qzRange'],
|
||||
qzRange=self._field_defaults['ReflectivityReductionConfig']['qzRange'],
|
||||
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):
|
||||
@@ -96,24 +96,24 @@ class FullAmorTest(TestCase):
|
||||
nu=0,
|
||||
muOffset=0.0,
|
||||
)
|
||||
reduction_config = options.ReductionConfig(
|
||||
normalisationMethod=self._field_defaults['ReductionConfig']['normalisationMethod'],
|
||||
reduction_config = options.ReflectivityReductionConfig(
|
||||
normalisationMethod=self._field_defaults['ReflectivityReductionConfig']['normalisationMethod'],
|
||||
qResolution=0.01,
|
||||
qzRange=self._field_defaults['ReductionConfig']['qzRange'],
|
||||
qzRange=self._field_defaults['ReflectivityReductionConfig']['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()
|
||||
|
||||
16
update.md
Normal file
16
update.md
Normal file
@@ -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`
|
||||
Reference in New Issue
Block a user