sinq-reflectometry/eos
0
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 5 Wiki Activity

implement events2histogram for LZ graph and automatic update

Browse Source
This commit is contained in:
glavic_a
2025-10-07 16:19:57 +02:00
parent 99af021b3e
commit 31201795b0
11 changed files with 235 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files
eos/e2h.py
+8 -4
View File
@@ -1,28 +1,32 @@
import logging
# need to do absolute import here as pyinstaller requires it
from eos.options import E2HConfig, ReaderConfig, ExperimentConfig
from eos.options import E2HConfig, ReaderConfig, ExperimentConfig, E2HReductionConfig
from eos.command_line import commandLineArgs
from eos.logconfig import setup_logging, update_loglevel
from eos.reduction_e2h import E2HReduction
def main():
setup_logging()
logging.getLogger('matplotlib').setLevel(logging.WARNING)
# read command line arguments and generate classes holding configuration parameters
clas = commandLineArgs([ReaderConfig, ExperimentConfig],
clas = commandLineArgs([ReaderConfig, ExperimentConfig, E2HReductionConfig],
'events2histogram')
update_loglevel(clas.verbose)
reader_config = ReaderConfig.from_args(clas)
experiment_config = ExperimentConfig.from_args(clas)
config = E2HConfig(reader_config, experiment_config, )
reduction_config = E2HReductionConfig.from_args(clas)
config = E2HConfig(reader_config, experiment_config, reduction_config)
logging.warning('######## events2histogram - data vizualization for Amor ########')
# only import heavy module if sufficient command line parameters were provided
from eos.reduction_reflectivity import ReflectivityReduction
# Create reducer with these arguments
reducer = ReflectivityReduction(config)
reducer = E2HReduction(config)
# Perform actual reduction
reducer.reduce()
eos/event_analysis.py
+1 -1
View File
@@ -18,7 +18,7 @@ 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
eos/event_data_types.py
+1 -1
View File
@@ -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)])
eos/event_handling.py
+3
View File
@@ -149,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
eos/file_reader.py
+15 -17
View File
@@ -30,7 +30,7 @@ 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"'
GREP = '/usr/bin/grep "value"'
else:
NICOS_CACHE_DIR = None
@@ -63,10 +63,11 @@ class AmorHeader:
try:
logging.warning(f" using parameter {nicos_key} from nicos cache")
year_date = self.fileDate.strftime('%Y')
value = str(subprocess.getoutput(f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year_date}')).split('\t')[-1]
call = f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year_date}'
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")
@@ -151,7 +152,7 @@ class AmorHeader:
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)
polarizationConfig = fileio.Polarization(polarizationConfigs[polarizationConfigLabel])
logging.debug(f' polarization configuration: {polarizationConfig} (index {polarizationConfigLabel})')
@@ -236,10 +237,6 @@ class AmorEventData(AmorHeader):
self.hdf = hdf
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()
@@ -253,12 +250,13 @@ class AmorEventData(AmorHeader):
"""
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]
@@ -279,15 +277,15 @@ class AmorEventData(AmorHeader):
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)
@@ -304,17 +302,17 @@ class AmorEventData(AmorHeader):
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):
eos/options.py
+8
View File
@@ -661,6 +661,14 @@ class E2HReductionConfig(ArgParsable):
},
)
plot_colormap: PlotColormaps = field(
default=PlotColormaps.gist_ncar,
metadata={
'short': 'pcmap',
'group': 'output',
'help': 'matplotlib colormap used in lambda-theta graphs when plotting',
},
)
@dataclass
eos/path_handling.py
+30 -4
View File
@@ -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,32 @@ 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}')
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()
return possible_indices[number-1]
eos/projection.py
+25 -2
View File
@@ -3,6 +3,8 @@ Classes used to calculate projections/binnings from event data onto given grids.
"""
import logging
from typing import Protocol
import numpy as np
from dataclasses import dataclass
@@ -75,6 +77,12 @@ class ProjectedReflectivity:
plt.xlabel('Q / $\\AA^{-1}$')
plt.ylabel('R')
class ProjectionInterface(Protocol):
def project(self, dataset: EventDatasetProtocol, monitor: float): ...
def clear(self): ...
def plot(self, **kwargs): ...
def update_plot(self): ...
class LZProjection:
grid: LZGrid
lamda: np.ndarray
@@ -165,6 +173,7 @@ class LZProjection:
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 +181,11 @@ 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
@property
def I(self):
output = self.data.I[:]
@@ -304,14 +318,23 @@ class LZProjection:
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)
self._graph = plt.pcolormesh(self.lamda, self.alphaF, self.data.ref, **kwargs)
if cmap:
plt.colorbar(label='R')
else:
if not 'norm' in kwargs:
kwargs['norm'] = LogNorm()
plt.pcolormesh(self.lamda, self.alphaF, self.data.I, **kwargs)
self._graph = plt.pcolormesh(self.lamda, self.alphaF, self.data.I, **kwargs)
if cmap:
plt.colorbar(label='I / cpm')
plt.xlabel('$\\lambda$ / $\\AA$')
plt.ylabel('$\\Theta$ / °')
def update_plot(self):
"""
Inline update of last plot by just updating the data.
"""
if self.is_normalized:
self._graph.set_array(self.data.ref)
else:
self._graph.set_array(self.data.I)
eos/reduction_e2h.py
+142 -11
View File
@@ -4,19 +4,29 @@ Can be used as a live preview with automatic update when files are modified.
"""
import logging
import os
import matplotlib.pyplot as plt
from .file_reader import AmorEventData
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, IncidentAngle, MonitorType, NormalisationMethod, MONITOR_UNITS, ExperimentConfig
from .options import E2HConfig, E2HPlotArguments, IncidentAngle, MonitorType, E2HPlotSelection
from . import event_handling as eh, event_analysis as ea
from .path_handling import PathResolver
from .projection import LZProjection, ProjectionInterface
class E2HReduction:
config: E2HConfig
header: Header
event_actions: eh.EventDataAction
_last_mtime = 0.
def __init__(self, config: E2HConfig):
self.config = config
@@ -28,23 +38,144 @@ class E2HReduction:
"""
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 = {}
self.fig = plt.figure()
if self.config.reduction.update:
# live update implies plotting
self.config.reduction.show_plot = True
# 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.ExtractWalltime()
if not self.config.reduction.fast:
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 or MonitorType.debug]:
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)
self.event_actions |= eh.FilterStrangeTimes()
self.event_actions |= ea.MergeFrames()
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)
if not self.config.reduction.fast:
self.event_actions |= eh.FilterStrangeTimes()
# select needed actions in depenence of plots
if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q,
E2HPlotSelection.L]:
self.event_actions |= ea.MergeFrames()
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()
self.grid = LZGrid([0.01], [0.0, 0.25])
# plot dependant options
if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]:
self.grid = LZGrid(0.01, [0.0, 0.25])
if self.config.reduction.plot==E2HPlotSelection.LT:
self.plot_kwds['colorbar'] = True
self.plot_kwds['cmap'] = str(self.config.reduction.plot_colormap)
def reduce(self):
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:
plt.savefig(f'e2h_{self.config.reduction.plot}.png', dpi=300)
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 prepare_graphs(self):
last_file_header = AmorHeader(self.file_list[-1])
tthh = last_file_header.geometry.nu - last_file_header.geometry.mu
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)
def read_data(self):
fileName = self.file_list[self.file_index]
self.file_index += 1
self.dataset = AmorEventData(fileName)
self.event_actions(self.dataset)
self.dataset.update_header(self.header)
self.header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1],
timestamp=self.dataset.fileDate))
def add_data(self):
self.monitor = self.dataset.data.pulses.monitor.sum()
self.projection.project(self.dataset, monitor=self.monitor)
def create_file_output(self):
...
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"\n$\\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 %I:%M:%S")
return output
def create_graph(self):
plt.title(self.create_title())
self.projection.plot(**self.plot_kwds)
plt.tight_layout()
def replace_dataset(self, latest):
self.file_list = self.path_resolver.resolve(f'{latest}')
self.file_index = 0
self.read_data()
self.projection.clear()
self.add_data()
self.fig.clear()
self.create_graph()
def update(self):
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 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)
except EOFError:
return
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.title(self.create_title())
plt.draw()
eos/reduction_reflectivity.py
+1 -1
View File
@@ -158,7 +158,7 @@ class ReflectivityReduction:
self.analyze_unsliced(i)
def analyze_unsliced(self, i):
self.monitor = np.sum(self.dataset.data.pulses.monitor)
self.monitor = self.dataset.data.pulses.monitor.sum()
logging.warning(f' monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.config.experiment.monitorType]}')
proj:LZProjection = self.project_on_lz()
setup.cfg
+1
View File
@@ -34,3 +34,4 @@ Homepage = "https://github.com/jochenstahn/amor"
[options.entry_points]
console_scripts =
eos = eos.__main__:main
events2histogram = eos.e2h:main