diff --git a/eos/normalization.py b/eos/normalization.py index 53cf570..9e64ed7 100644 --- a/eos/normalization.py +++ b/eos/normalization.py @@ -66,6 +66,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) diff --git a/eos/options.py b/eos/options.py index 5d06d14..3b798f0 100644 --- a/eos/options.py +++ b/eos/options.py @@ -662,6 +662,15 @@ class E2HReductionConfig(ArgParsable): }, ) + 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.gist_ncar, metadata={ diff --git a/eos/projection.py b/eos/projection.py index fbe0b74..9b2753d 100644 --- a/eos/projection.py +++ b/eos/projection.py @@ -4,6 +4,7 @@ Classes used to calculate projections/binnings from event data onto given grids. import logging from abc import ABC, abstractmethod +from typing import List, Tuple, Union import numpy as np from dataclasses import dataclass @@ -26,7 +27,7 @@ class ProjectionInterface(ABC): def update_plot(self): ... @dataclass -class ProjectedReflectivity(ProjectionInterface): +class ProjectedReflectivity: R: np.ndarray dR: np.ndarray Q: np.ndarray @@ -83,33 +84,6 @@ class ProjectedReflectivity(ProjectionInterface): 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 - 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') - - def update_plot(self): - ln, (errx_top, errx_bot, erry_top, erry_bot), (barsx, barsy) = self._graph.lines - - 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 @@ -318,6 +292,7 @@ class LZProjection(ProjectionInterface): plt.xlim(3., 12.) 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) @@ -339,14 +314,57 @@ class LZProjection(ProjectionInterface): 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') + 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(plt.gca().lines)+list(plt.gca().texts): + 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.yscale('log') + plt.xlabel('Q / $\\AA^{-1}$') + plt.ylabel('R') + + def update_plot(self): + ln, (errx_top, errx_bot, erry_top, erry_bot), (barsx, barsy) = self._graph.lines + + yerr_top = self.data.R+self.data.dR + yerr_bot = self.data.R-self.data.dR + + errx_top.set_ydata(self.data.R) + errx_bot.set_ydata(self.data.R) + erry_top.set_ydata(yerr_top) + erry_bot.set_ydata(yerr_bot) + + ln.set_ydata(self.data.R) + class YZProjection(ProjectionInterface): y: np.ndarray @@ -405,6 +423,7 @@ class YZProjection(ProjectionInterface): plt.ylabel('Z') plt.xlim(self.y[0], self.y[-1]) plt.ylim(self.z[0], self.z[-1]) + plt.title('Horizontal Pixel vs. Vertical Pixel') self._graph_axis = plt.gca() plt.connect('button_press_event', self.draw_yzcross) @@ -421,12 +440,12 @@ class YZProjection(ProjectionInterface): 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') + 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(plt.gca().lines)+list(plt.gca().texts): + for art in list(self._graph_axis.lines)+list(self._graph_axis.texts): art.remove() plt.draw() @@ -490,6 +509,7 @@ class TofZProjection(ProjectionInterface): plt.ylabel('Z') plt.xlim(self.tof[0]*1e3, self.tof[-1]*1e3) plt.ylim(self.z[0], self.z[-1]) + plt.title('Time of Flight vs. Vertical Pixel') self._graph_axis = plt.gca() plt.connect('button_press_event', self.draw_tzcross) @@ -506,11 +526,49 @@ class TofZProjection(ProjectionInterface): 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.tof[0]*1e3, self.tof[-1]*1e3], [event.ydata, event.ydata], '-', color='grey') - plt.text(event.xdata, event.ydata, f'({event.xdata:.2f}, {event.ydata:.1f})', backgroundcolor='white') + 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(plt.gca().lines)+list(plt.gca().texts): + for art in list(self._graph_axis.lines)+list(self._graph_axis.texts): art.remove() plt.draw() + +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]) + 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)) + + def update_plot(self): + for pi in self.projections: + pi.update_plot() diff --git a/eos/reduction_e2h.py b/eos/reduction_e2h.py index 75b8c76..7aaab28 100644 --- a/eos/reduction_e2h.py +++ b/eos/reduction_e2h.py @@ -18,7 +18,8 @@ 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 LZProjection, ProjectionInterface, TofZProjection, YZProjection +from .projection import CombinedProjection, LZProjection, ProjectionInterface, ReflectivityProjector, TofZProjection, \ + YZProjection NEEDS_LAMDA = (E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q, E2HPlotSelection.L) @@ -85,13 +86,16 @@ 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 in [E2HPlotSelection.LT, E2HPlotSelection.YZ, E2HPlotSelection.TZ]: + if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.YZ, E2HPlotSelection.TZ]: self.plot_kwds['colorbar'] = True self.plot_kwds['cmap'] = str(self.config.reduction.plot_colormap) def reduce(self): if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]: - self.norm = LZNormalisation.unity(self.grid) + if self.config.reduction.normalizationModel: + self.norm = LZNormalisation.model(self.grid) + else: + self.norm = LZNormalisation.unity(self.grid) self.prepare_graphs() @@ -123,6 +127,17 @@ class E2HReduction: 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_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_norm_mask(self.norm) + self.projection = ReflectivityProjector(plz, self.norm) if self.config.reduction.plot==E2HPlotSelection.YZ: self.projection = YZProjection() @@ -130,6 +145,18 @@ class E2HReduction: 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.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_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)]) + def read_data(self): fileName = self.file_list[self.file_index] self.file_index += 1 @@ -143,6 +170,8 @@ class E2HReduction: 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): ... @@ -157,7 +186,7 @@ class E2HReduction: return output def create_graph(self): - plt.title(self.create_title()) + plt.suptitle(self.create_title()) self.projection.plot(**self.plot_kwds) plt.tight_layout()