From 3e6c121aeb2db4419841387910a04041bb7f426e Mon Sep 17 00:00:00 2001 From: Philipp Niedermayer <p.niedermayer@gsi.de> Date: Fri, 19 Aug 2022 14:51:55 +0200 Subject: [PATCH] Restructure module --- plotting/__init__.py | 18 ++++ plotting/common.py | 60 ++++++++++++ plotting.py => plotting/other.py | 155 +------------------------------ plotting/spill.py | 140 ++++++++++++++++++++++++++++ 4 files changed, 219 insertions(+), 154 deletions(-) create mode 100644 plotting/__init__.py create mode 100644 plotting/common.py rename plotting.py => plotting/other.py (83%) create mode 100644 plotting/spill.py diff --git a/plotting/__init__.py b/plotting/__init__.py new file mode 100644 index 0000000..b86897b --- /dev/null +++ b/plotting/__init__.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" Methos for plotting measurement data + +""" + +__author__ = "Philipp Niedermayer" +__contact__ = "p.niedermayer@gsi.de" +__date__ = "2022-05-04" + + +from .spill import * +from .other import * + + + + diff --git a/plotting/common.py b/plotting/common.py new file mode 100644 index 0000000..0f4212f --- /dev/null +++ b/plotting/common.py @@ -0,0 +1,60 @@ +import scipy +import scipy.signal +import numpy as np +import matplotlib as mpl +import matplotlib.pyplot as plt +import pint + +from ..data import * +from ..fitting import * + + +# Setup pint +SI = pint.UnitRegistry() +SI.setup_matplotlib() +SI.default_format = '~P' + + +# Colors +# https://arxiv.org/abs/2107.02270 +petroff_colors = ["#3f90da", "#ffa90e", "#bd1f01", "#94a4a2", "#832db6", "#a96b59", "#e76300", "#b9ac70", "#717581", "#92dadd"] +cmap_petroff_10 = mpl.colors.ListedColormap(petroff_colors, 'Petroff 10') +mpl.rcParams['axes.prop_cycle'] = mpl.cycler(color=petroff_colors) + +cmap_petroff_gradient = mpl.colors.LinearSegmentedColormap.from_list('Petroff gradient', [petroff_colors[i] for i in (9,0,4,2,6,1)]) +cmap_petroff_gradient.set_under(petroff_colors[3]) +cmap_petroff_gradient.set_over(petroff_colors[7]) +mpl.rcParams['image.cmap'] = cmap_petroff_gradient + +cmap_petroff_bipolar = mpl.colors.LinearSegmentedColormap.from_list('Petroff bipolar', [petroff_colors[i] for i in (2,6,1,3,9,0,4)]) +cmap_petroff_bipolar.set_under(petroff_colors[5]) +cmap_petroff_bipolar.set_over(petroff_colors[8]) + + +# Matplotlib options +mpl.rcParams.update({ + 'figure.figsize': (8, 5), # Note: overwritten when using the %matplotlib magics + 'figure.dpi': 120, # Note: overwritten when using the %matplotlib magics + 'figure.constrained_layout.use': True, + 'legend.fontsize': 'x-small', + 'legend.title_fontsize': 'small', + 'grid.color': '#DDD', +}) + + + + +# Utility methods + + + +def ax_set(ax, swap_xy=False, **kwargs): + if swap_xy: + nkwargs = {} + for k in kwargs: + nk = (('y' if k[0]=='x' else 'x') + k[1:]) if k[0] in 'xy' else k + nkwargs[nk] = kwargs[k] + ax.set(**nkwargs) + else: + ax.set(**kwargs) + diff --git a/plotting.py b/plotting/other.py similarity index 83% rename from plotting.py rename to plotting/other.py index 9890cb5..695cb3f 100644 --- a/plotting.py +++ b/plotting/other.py @@ -1,46 +1,7 @@ -import scipy -import scipy.signal -import numpy as np -import matplotlib as mpl -import matplotlib.pyplot as plt -import pint -from .data import * -from .fitting import * +from .common import * -# Setup pint -SI = pint.UnitRegistry() -SI.setup_matplotlib() -SI.default_format = '~P' - - -# Colors -# https://arxiv.org/abs/2107.02270 -petroff_colors = ["#3f90da", "#ffa90e", "#bd1f01", "#94a4a2", "#832db6", "#a96b59", "#e76300", "#b9ac70", "#717581", "#92dadd"] -cmap_petroff_10 = mpl.colors.ListedColormap(petroff_colors, 'Petroff 10') -mpl.rcParams['axes.prop_cycle'] = mpl.cycler(color=petroff_colors) - -cmap_petroff_gradient = mpl.colors.LinearSegmentedColormap.from_list('Petroff gradient', [petroff_colors[i] for i in (9,0,4,2,6,1)]) -cmap_petroff_gradient.set_under(petroff_colors[3]) -cmap_petroff_gradient.set_over(petroff_colors[7]) -mpl.rcParams['image.cmap'] = cmap_petroff_gradient - -cmap_petroff_bipolar = mpl.colors.LinearSegmentedColormap.from_list('Petroff bipolar', [petroff_colors[i] for i in (2,6,1,3,9,0,4)]) -cmap_petroff_bipolar.set_under(petroff_colors[5]) -cmap_petroff_bipolar.set_over(petroff_colors[8]) - - -# Matplotlib options -mpl.rcParams.update({ - 'figure.figsize': (8, 5), # Note: overwritten when using the %matplotlib magics - 'figure.dpi': 120, # Note: overwritten when using the %matplotlib magics - 'figure.constrained_layout.use': True, - 'legend.fontsize': 'x-small', - 'legend.title_fontsize': 'small', - 'grid.color': '#DDD', -}) - def add_vline(axes, r, color='k', text=None, label=None, order=None, lw=1, text_top=True, text_vertical=True, zorder=-100, **kwargs): @@ -240,15 +201,6 @@ def add_scale(ax, scale, text=None, *, vertical=False, size=0.01, padding=0.1, l def v(swap_xy, x, y): return (y, x) if swap_xy else (x, y) -def ax_set(ax, swap_xy=False, **kwargs): - if swap_xy: - nkwargs = {} - for k in kwargs: - nk = (('y' if k[0]=='x' else 'x') + k[1:]) if k[0] in 'xy' else k - nkwargs[nk] = kwargs[k] - ax.set(**nkwargs) - else: - ax.set(**kwargs) def smooth_plot(ax, x, y, smoothing=None, swap_xy=False, **kwargs): @@ -760,108 +712,3 @@ def plot_beam_size(ax, ipm_data, xy, time_range=None, smoothing=None, **kwargs): -# TDC data -########### - -def plot_spill_intensity(ax, spill_data, bin_time=10e-6, *, rate=False, cumulative=False, relative=False, **plot_kwargs): - """Plot spill intensity (binned particle counts) over extraction time - - :param ax: Axis to plot onto - :param spill_data: Instance of TDCSpill - :param bin_time: width of bins for particle counting in sec - :param rate: if true, plot particle count rate, i.e. counts per s instead of per bin - :param cumulative: if true, plot cumulative particle number - :param relative: if true, plot fraction of total particle number - """ - nbins = int(np.ceil(spill_data.length/bin_time)) - - assert not (rate and cumulative), 'rate and cumulative plotting are mutually exclusive' - weights = np.ones_like(spill_data.hittimes) - if rate: weights /= bin_time - if relative: weights /= spill_data.counts - label = 'Extracted particle' + \ - (' fraction' if relative else 's') + \ - ('' if cumulative else \ - (' / s' if rate else f'\nper ${(bin_time*SI.s).to_compact():~gL}$ interval')) - - ax.hist(spill_data.hittimes, range=(spill_data.time_offset, spill_data.time_offset+bin_time*nbins), bins=nbins, - histtype='step', cumulative=cumulative, weights=weights, **plot_kwargs) - ax.set(ylabel=label, xlabel='Extraction time / s', #xlim=(spill_data.time_offset, spill_data.time_offset+bin_time*nbins), - ) - if relative: ax.yaxis.set_major_formatter(mpl.ticker.PercentFormatter(1)) - - -def plot_spill_duty_factor(ax, spill_data, counting_dt=10e-6, evaluation_dt=10e-3, *, show_poisson_limit=False, **kwargs): - """Plot spill duty factor over extraction time - - :param ax: Axis to plot onto - :param spill_data: Instance of TDCSpill - :param counting_dt: width of bins for particle counting in sec - :param evaluation_dt: width of bins for duty factor evaluation (value is rounded to the next multiple of count_bin_time) - - :param show_poisson_limit: if true, show poison limit of duty factor - - """ - - counts, edges = spill_data.hitcounts(counting_dt, evaluation_dt) - - # spill duty factor - F = np.mean(counts, axis=1)**2 / np.mean(counts**2, axis=1) - s = ax.stairs(F, edges, **kwargs) - - if show_poisson_limit: - Fp = np.mean(counts, axis=1) / ( np.mean(counts, axis=1) + 1 ) - ax.stairs(Fp, edges, label='Poisson limit', color=s.get_edgecolor() or 'gray', alpha=0.5,zorder=-1, lw=1, ls=':') - - ax.set(xlabel='Extraction time / s', ylim=(0,1.1), - ylabel=f'Spill duty factor F\n(${(evaluation_dt*SI.s).to_compact():~gL}$ / ${(counting_dt*SI.s).to_compact():~gL}$ intervals)' - ) - - -def plot_spill_consecutive_particle_delay(ax, spill_data, limit=5e-6, resolution=10e-9, show_poisson=False, **plot_kwargs): - """Plot consecutive particle separation (particle-interval histograms) - - :param ax: Axis to plot onto - :param spill_data: Instance of TDCSpill - :param limit: maximum delay to plot (range of histogram) - :param resolution: time resolution (bin width of histogram) - """ - nbins = int(np.ceil(limit/resolution)) - delay = spill_data.hitdelay - hist = _, t = np.histogram(delay, range=(0, resolution*nbins), bins=nbins) #, weights=np.ones_like(delay)/resolution) - ax.stairs(*hist, **plot_kwargs) - ax.set(xlabel='Delay between consecutive particles / s', xlim=(resolution, resolution*nbins), - ylabel=f'particle count / ${(resolution*SI.s).to_compact():~gL}$ delay interval') - - if show_poisson: - # plot poisson statistics (exponential distribution) - l = spill_data.counts/spill_data.length - p = np.exp(-l*t[:-1]) - np.exp(-l*t[1:]) # probability of delay in between t[i] and t[i+1] - c = p*spill_data.counts - c[c<1] = 0 # supress fractional counts - ax.stairs(c, t, color='gray', label='Poisson statistics', zorder=-1, lw=1) - - -def plot_spill_frequency_spectrum(ax, spill_data, fmax=100e3, **plot_kwargs): - """Plot frequency spectrum of particle arrival times - - :param ax: Axis to plot onto - :param spill_data: Instance of TDCSpill - :param fmax: maximum frequency for analysis in Hz - :param bins: number of bins for fft - """ - - # re-sample timestamps into equally sampled time series - dt = 1/2/fmax - timeseries, _ = spill_data.hitcounts(dt) - - # frequency analysis - freq = np.fft.rfftfreq(len(timeseries), d=dt)[1:] - mag = np.abs(np.fft.rfft(timeseries))[1:] - - # plot - ax.plot(freq, mag, **plot_kwargs) - ax.set(xlabel=f'Frequency / Hz', xlim=(0, fmax), - ylabel='Magnitude / a.u.', ylim=(0, 1.1*max(mag))) - - diff --git a/plotting/spill.py b/plotting/spill.py new file mode 100644 index 0000000..178d960 --- /dev/null +++ b/plotting/spill.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" Methos for making spill realted plots + +""" + +__author__ = "Philipp Niedermayer" +__contact__ = "p.niedermayer@gsi.de" +__date__ = "2022-08-19" + + + +from .common import * + + + + +def plot_spill_intensity(ax, spill_data, bin_time=10e-6, *, rate=False, cumulative=False, relative=False, **plot_kwargs): + """Plot spill intensity (binned particle counts) over extraction time + + :param ax: Axis to plot onto + :param spill_data: Instance of spill data class, e.g. TDCSpill + :param bin_time: Width of bins for particle counting in sec + :param rate: If true, plot particle count rate, i.e. counts per s instead of per bin + :param cumulative: If true, plot cumulative particle number + :param relative: If true, plot fraction of total particle number + :param plot_kwargs: Keyword arguments to be passed to plotting method + """ + assert not (rate and cumulative), 'rate and cumulative plotting are mutually exclusive' + + # histogram in time bins + nbins = int(np.ceil(spill_data.length/bin_time)) + weights = np.ones_like(spill_data.hittimes) + if rate: weights /= bin_time + if relative: weights /= spill_data.counts + hist, edges = np.histogram(spill_data.hittimes, bins=nbins, weights=weights, + range=(spill_data.time_offset, spill_data.time_offset+bin_time*nbins)) + if cumulative: hist = np.cumsum(hist) + + # plot + ax.stairs(hist, edges, **plot_kwargs) + if relative: ax.yaxis.set_major_formatter(mpl.ticker.PercentFormatter(1)) + ax.set(xlabel='Extraction time / s', + ylabel='Extracted particle' + \ + (' fraction' if relative else 's') + \ + ('' if cumulative else \ + (' / s' if rate else f'\nper ${(bin_time*SI.s).to_compact():~gL}$ interval')), + ) + + + + +def plot_spill_duty_factor(ax, spill_data, counting_dt=10e-6, evaluation_dt=10e-3, *, show_poisson_limit=False, **plot_kwargs): + """Plot spill duty factor (<x>²/<x²>) over extraction time + + :param ax: Axis to plot onto + :param spill_data: Instance of spill data class, e.g. TDCSpill + :param counting_dt: Width of bins for particle counting in sec + :param evaluation_dt: Width of bins for duty factor evaluation (value is rounded to the next multiple of count_bin_time) + :param show_poisson_limit: If true, show poison limit of duty factor + :param plot_kwargs: Keyword arguments to be passed to plotting method + """ + + counts, edges = spill_data.hitcounts(counting_dt, evaluation_dt) + + # spill duty factor + F = np.mean(counts, axis=1)**2 / np.mean(counts**2, axis=1) + s = ax.stairs(F, edges, **plot_kwargs) + + if show_poisson_limit: + Fp = np.mean(counts, axis=1) / ( np.mean(counts, axis=1) + 1 ) + ax.stairs(Fp, edges, label='Poisson limit', color=s.get_edgecolor() or 'gray', alpha=0.5,zorder=-1, lw=1, ls=':') + + ax.set(xlabel='Extraction time / s', ylim=(0,1.1), + ylabel=f'Spill duty factor F\n(${(evaluation_dt*SI.s).to_compact():~gL}$ / ${(counting_dt*SI.s).to_compact():~gL}$ intervals)' + ) + + + + +def plot_spill_consecutive_particle_delay(ax, spill_data, limit=5e-6, resolution=10e-9, show_poisson=False, swap_axis=False, log=True, + **plot_kwargs): + """Plot histogram of consecutive particle separation (particle-interval histograms) + + :param ax: Axis to plot onto + :param spill_data: Instance of TDCSpill + :param limit: Maximum delay to plot (range of histogram) + :param resolution: Time resolution (bin width of histogram) + :param show_poisson: If true, show poison like distribution + :param log: If true, use logarithmic scaling + :param swap_axis: If True, swap plot axis + :param plot_kwargs: Keyword arguments to be passed to plotting method + """ + nbins = int(np.ceil(limit/resolution)) + delay = spill_data.hitdelay + hist = c, t = np.histogram(delay, range=(0, resolution*nbins), bins=nbins) #, weights=np.ones_like(delay)/resolution) + ax.stairs(*hist, orientation='horizontal' if swap_axis else 'vertical', **plot_kwargs) + ax_set(ax, swap_axis, + xlabel='Delay between consecutive particles / s', xlim=(resolution, resolution*nbins), xscale='log' if log else None, + ylabel=f'Particle count / ${(resolution*SI.s).to_compact():~gL}$ delay interval', yscale='log' if log else None) + + if show_poisson: + # plot poisson statistics (exponential distribution) + l = spill_data.counts/spill_data.length + p = np.exp(-l*t[:-1]) - np.exp(-l*t[1:]) # probability of delay in between t[i] and t[i+1] + c = p*spill_data.counts + c[c<1] = 0 # supress fractional counts + ax.stairs(c, t, orientation='horizontal' if swap_axis else 'vertical', color='gray', label='Poisson statistics', zorder=-1, lw=1) + + + + +def plot_spill_frequency_spectrum(ax, spill_data, fmax=100e3, **plot_kwargs): + """Plot frequency spectrum of particle arrival times + + :param ax: Axis to plot onto + :param spill_data: Instance of TDCSpill + :param fmax: Maximum frequency for analysis in Hz + :param plot_kwargs: Keyword arguments to be passed to plotting method + """ + + # re-sample timestamps into equally sampled time series + dt = 1/2/fmax + timeseries, _ = spill_data.hitcounts(dt) + + # frequency analysis + freq = np.fft.rfftfreq(len(timeseries), d=dt)[1:] + mag = np.abs(np.fft.rfft(timeseries))[1:] + + # plot + ax.plot(freq, mag, **plot_kwargs) + ax.set(xlabel=f'Frequency / Hz', xlim=(0, fmax), + ylabel='Magnitude / a.u.', ylim=(0, 1.1*max(mag))) + + + + + + \ No newline at end of file -- GitLab