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

Add YZProjection

Browse Source
This commit is contained in:
glavic_a
2025-10-08 08:54:39 +02:00
parent 2972ffd5d8
commit b71b72c6a3
3 changed files with 170 additions and 84 deletions
Download Patch File Download Diff File Expand all files Collapse all files
eos/instrument.py
+1
View File
@@ -99,6 +99,7 @@ class LZGrid:
@cache
def z(self):
# TODO: shouldn't this be -0.5 to be the edges of each pixel?
return np.arange(Detector.nBlades*Detector.nWires+1)
@cache
eos/projection.py
+152 -64
View File
@@ -3,7 +3,7 @@ Classes used to calculate projections/binnings from event data onto given grids.
"""
import logging
from typing import Protocol
from abc import ABC, abstractmethod
import numpy as np
from dataclasses import dataclass
@@ -12,8 +12,21 @@ 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:
class ProjectedReflectivity(ProjectionInterface):
R: np.ndarray
dR: np.ndarray
Q: np.ndarray
@@ -70,26 +83,51 @@ class ProjectedReflectivity:
self.R -= R
self.dR = np.sqrt(self.dR**2+dR**2)
def project(self, dataset: EventDatasetProtocol, monitor: float):
raise NotImplementedError("Direct projection from dataset is not yet available")
def clear(self):
raise NotImplementedError("Direct projection from dataset is not yet available")
def plot(self, **kwargs):
from matplotlib import pyplot as plt
plt.errorbar(self.Q, self.R, xerr=self.dQ, yerr=self.dR, **kwargs)
self._graph = plt.errorbar(self.Q, self.R, xerr=self.dQ, yerr=self.dR, **kwargs)
self._graph_axis = plt.gca()
plt.yscale('log')
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): ...
def update_plot(self):
ln, (errx_top, errx_bot, erry_top, erry_bot), (barsx, barsy) = self._graph.lines
class LZProjection:
yerr_top = self.R+self.dR
yerr_bot = self.R-self.dR
errx_top.set_ydata(self.R)
errx_bot.set_ydata(self.R)
erry_top.set_ydata(yerr_top)
erry_bot.set_ydata(yerr_bot)
ln.set_ydata(self.R)
class LZProjection(ProjectionInterface):
grid: LZGrid
lamda: np.ndarray
alphaF: np.ndarray
is_normalized: bool
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
@@ -103,16 +141,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.
@@ -185,6 +214,7 @@ class LZProjection:
# empty data
self.data[:] = 0
self.data.mask = True
self.monitor = 0.
@property
def I(self):
@@ -261,50 +291,6 @@ class LZProjection:
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
@@ -329,12 +315,114 @@ class LZProjection:
plt.colorbar(label='I / cpm')
plt.xlabel('$\\lambda$ / $\\AA$')
plt.ylabel('$\\Theta$ / °')
plt.xlim(3., 12.)
af = self.alphaF[self.data.mask]
plt.ylim(af.min(), af.max())
self._graph_axis = plt.gca()
plt.connect('button_press_event', self.draw_qline)
def update_plot(self):
"""
Inline update of last plot by just updating the data.
Inline update of previous 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)
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
plt.plot([0, xmax], [0, slope*xmax], '-', color='grey')
plt.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(plt.gca().lines)+list(plt.gca().texts):
art.remove()
plt.draw()
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
if not 'aspect' in kwargs:
kwargs['aspect'] = 'auto'
if not 'norm' in kwargs:
kwargs['norm'] = LogNorm()
self._graph = plt.imshow(self.data.I[:, ::-1].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[0], self.z[-1])
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.
"""
self._graph.set_array(self.data.I[:, ::-1].T)
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=='':
plt.plot([event.xdata, event.xdata], [self.z[0], self.z[-1]], '-', color='grey')
plt.plot([self.y[0], self.y[-1]], [event.ydata, event.ydata], '-', color='grey')
plt.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(plt.gca().lines)+list(plt.gca().texts):
art.remove()
plt.draw()
eos/reduction_e2h.py
+17 -20
View File
@@ -18,7 +18,7 @@ from .normalization import LZNormalisation
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
from .projection import LZProjection, ProjectionInterface, YZProjection
class E2HReduction:
@@ -27,6 +27,7 @@ class E2HReduction:
event_actions: eh.EventDataAction
_last_mtime = 0.
projection: ProjectionInterface
def __init__(self, config: E2HConfig):
self.config = config
@@ -76,12 +77,13 @@ class E2HReduction:
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:
if self.config.reduction.plot in [E2HPlotSelection.LT, E2HPlotSelection.YZ]:
self.plot_kwds['colorbar'] = True
self.plot_kwds['cmap'] = str(self.config.reduction.plot_colormap)
def reduce(self):
self.norm = LZNormalisation.unity(self.grid)
if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]:
self.norm = LZNormalisation.unity(self.grid)
self.prepare_graphs()
@@ -94,7 +96,8 @@ class E2HReduction:
if self.config.reduction.plotArgs!=E2HPlotArguments.OutputFile or self.config.reduction.show_plot:
self.create_graph()
if self.config.reduction.plotArgs==E2HPlotArguments.Default:
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.update:
self.timer = self.fig.canvas.new_timer(1000)
@@ -113,6 +116,9 @@ class E2HReduction:
if not self.config.reduction.fast:
self.projection.correct_gravity(last_file_header.geometry.detectorDistance)
if self.config.reduction.plot==E2HPlotSelection.YZ:
self.projection = YZProjection()
def read_data(self):
fileName = self.file_list[self.file_index]
self.file_index += 1
@@ -143,23 +149,13 @@ class E2HReduction:
plt.title(self.create_title())
self.projection.plot(**self.plot_kwds)
plt.tight_layout()
if self.config.reduction.plot==E2HPlotSelection.LT:
plt.connect('button_press_event', self.draw_qline)
def draw_qline(self, event):
if event.button is plt.MouseButton.LEFT and self.fig.canvas.manager.toolbar.mode=='':
slope = event.ydata/event.xdata
xmax = 12.5
plt.plot([0, xmax], [0, slope*xmax], '-', color='grey')
plt.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 self.fig.canvas.manager.toolbar.mode=='':
for art in list(plt.gca().lines)+list(plt.gca().texts):
art.remove()
plt.draw()
def replace_dataset(self, latest):
self.file_list = self.path_resolver.resolve(f'{latest}')
new_files = self.path_resolver.resolve(f'{latest}')
if not os.path.exists(new_files[-1]):
return
logging.warning(f"Preceding to next file {latest}")
self.file_list = new_files
self.file_index = 0
self.read_data()
self.projection.clear()
@@ -168,6 +164,7 @@ class E2HReduction:
self.create_graph()
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])
@@ -183,7 +180,7 @@ class E2HReduction:
update_data = AmorEventData(self.file_list[-1], self.dataset.last_index+1)
except EOFError:
return
logging.info(" updating with new data")
self.event_actions(update_data)
self.dataset=update_data