In [1]:
import os, math
from dataclasses import dataclass
import numpy as np
import pandas as pd
from numba import vectorize, float64, boolean, njit
from scipy.signal import find_peaks
import matplotlib.pyplot as plt
import holoviews as hv
from IPython.display import Markdown
from holoviews.operation import decimate
from holoviews.operation.datashader import datashade,rasterize, dynspread
import dask.dataframe as dd
# import hvplot.pandas
# import hvplot.dask
hv.extension('bokeh', logo = False)
In [2]:
# turn off numpy warning for overflow in np.exp()
np.seterr(over='ignore');
In [3]:
shot_no = 0
(Internal Tektronix format) .wfm parser¶
In [4]:
# wfm reader proof-of-concept
# https://www.tek.com/sample-license
# reads volts vs. time records (including fastframes) from little-endian version 3 WFM files
# See Also
# Performance Oscilloscope Reference Waveform File Format
# Tektronix part # 077-0220-10
# https://www.tek.com/oscilloscope/dpo7000-digital-phosphor-oscilloscope-manual-4
import struct
import time
import numpy as np # http://www.numpy.org/
class WfmReadError(Exception):
"""error for unexpected things"""
pass
def read_wfm(target):
"""return sample data from target WFM file"""
with open(target, 'rb') as f:
hbytes = f.read(838)
meta = decode_header(hbytes)
# file signature checks
if meta['byte_order'] != 0x0f0f:
raise WfmReadError('big-endian not supported in this example')
if meta['version'] != b':WFM#003':
raise WfmReadError('only version 3 wfms supported in this example')
if meta['imp_dim_count'] != 1:
raise WfmReadError('imp dim count not 1')
if meta['exp_dim_count'] != 1:
raise WfmReadError('exp dim count not 1')
if meta['record_type'] != 2:
raise WfmReadError('not WFMDATA_VECTOR')
if meta['exp_dim_1_type'] != 0:
raise WfmReadError('not EXPLICIT_SAMPLE')
if meta['time_base_1'] != 0:
raise WfmReadError('not BASE_TIME')
tfrac_array = np.zeros(meta['Frames'], dtype=np.double)
tdatefrac_array = np.zeros(meta['Frames'], dtype=np.double)
tdate_array = np.zeros(meta['Frames'], dtype=np.int32)
tfrac_array[0] = meta['tfrac']
tdatefrac_array[0] = meta['tdatefrac']
tdate_array[0] = meta['tdate']
# if fastframe, read fastframe table
if meta['fastframe'] == 1:
WUSp = np.fromfile(f, dtype='i4,f8,f8,i4', count=(meta['Frames'] - 1))
# merge first frame trigger infos with frames > 1
tfrac_array[1:] = WUSp['f1']
tdatefrac_array[1:] = WUSp['f2']
tdate_array[1:] = WUSp['f3']
# read curve block
bin_wave = np.memmap(filename = f,
dtype = meta['dformat'],
mode = 'r',
offset = meta['curve_offset'],
shape = (meta['avilable_values'], meta['Frames']),
order = 'F')
# close file
# slice out buffer values
bin_wave = bin_wave[meta['pre_values']:meta['avilable_values'] - meta['post_values'],:]
scaled_array = bin_wave * meta['vscale'] + meta['voffset']
return scaled_array, meta['tstart'], meta['tscale'], tfrac_array, tdatefrac_array, tdate_array
def decode_header(header_bytes):
"""returns a dict of wfm metadata"""
wfm_info = {}
if len(header_bytes) != 838:
raise WfmReadError('wfm header bytes not 838')
wfm_info['byte_order'] = struct.unpack_from('H', header_bytes, offset=0)[0]
wfm_info['version'] = struct.unpack_from('8s', header_bytes, offset=2)[0]
wfm_info['imp_dim_count'] = struct.unpack_from('I', header_bytes, offset=114)[0]
wfm_info['exp_dim_count'] = struct.unpack_from('I', header_bytes, offset=118)[0]
wfm_info['record_type'] = struct.unpack_from('I', header_bytes, offset=122)[0]
wfm_info['exp_dim_1_type'] = struct.unpack_from('I', header_bytes, offset=244)[0]
wfm_info['time_base_1'] = struct.unpack_from('I', header_bytes, offset=768)[0]
wfm_info['fastframe'] = struct.unpack_from('I', header_bytes, offset=78)[0]
wfm_info['Frames'] = struct.unpack_from('I', header_bytes, offset=72)[0] + 1
wfm_info['summary_frame'] = struct.unpack_from('h', header_bytes, offset=154)[0]
wfm_info['curve_offset'] = struct.unpack_from('i', header_bytes, offset=16)[0] # 838 + ((frames - 1) * 54)
# scaling factors
wfm_info['vscale'] = struct.unpack_from('d', header_bytes, offset=168)[0]
wfm_info['voffset'] = struct.unpack_from('d', header_bytes, offset=176)[0]
wfm_info['tstart'] = struct.unpack_from('d', header_bytes, offset=496)[0]
wfm_info['tscale'] = struct.unpack_from('d', header_bytes, offset=488)[0]
# trigger detail
wfm_info['tfrac'] = struct.unpack_from('d', header_bytes, offset=788)[0] # frame index 0
wfm_info['tdatefrac'] = struct.unpack_from('d', header_bytes, offset=796)[0] # frame index 0
wfm_info['tdate'] = struct.unpack_from('I', header_bytes, offset=804)[0] # frame index 0
# data offsets
# frames are same size, only first frame offsets are used
dpre = struct.unpack_from('I', header_bytes, offset=822)[0]
wfm_info['dpre'] = dpre
dpost = struct.unpack_from('I', header_bytes, offset=826)[0]
wfm_info['dpost'] = dpost
readbytes = dpost - dpre
wfm_info['readbytes'] = readbytes
allbytes = struct.unpack_from('I', header_bytes, offset=830)[0]
wfm_info['allbytes'] = allbytes
# sample data type detection
code = struct.unpack_from('i', header_bytes, offset=240)[0]
wfm_info['code'] = code
bps = struct.unpack_from('b', header_bytes, offset=15)[0] # bytes-per-sample
wfm_info['bps'] = bps
if code == 7 and bps == 1:
dformat = 'int8'
samples = readbytes
elif code == 0 and bps == 2:
dformat = 'int16'
samples = readbytes // 2
elif code == 4 and bps == 4:
dformat = 'single'
samples = readbytes // 4
else:
raise WfmReadError('data type code or bytes-per-sample not understood')
wfm_info['dformat'] = dformat
wfm_info['samples'] = samples
wfm_info['avilable_values'] = allbytes // bps
wfm_info['pre_values'] = dpre // bps
wfm_info['post_values'] = (allbytes - dpost) // bps
return wfm_info
def readwfm(path):
volts, tstart, tscale, tfrac, tdatefrac, tdate = read_wfm(path)
toff = tfrac * tscale
samples, frames = volts.shape
tstop = samples * tscale + tstart
volts = volts.reshape(len(volts))
time = np.linspace(tstart+toff, tstop+toff, num=samples, endpoint=False)
time = time.reshape(len(volts))
return time,volts
Get Hi-Res data from MSO58 oscilloscope¶
In [5]:
# load data locally
fnames = ['DAS_raw_data_dir/ch2.wfm',
'DAS_raw_data_dir/ch3.wfm',
'DAS_raw_data_dir/ch4.wfm',
'DAS_raw_data_dir/ch5.wfm',
'DAS_raw_data_dir/ch6.wfm',
'DAS_raw_data_dir/ch7.wfm']
ds = np.DataSource('/tmp/')
if not os.path.isfile('DAS_raw_data_dir/ch6.wfm'):
# or download from web (slow)
for i, ch_no in enumerate(range(2,8)):
try:
file = ds.open(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Devices/Oscilloscopes/TektrMSO58-a/ch{ch_no}.wfm', mode = 'rb')
fnames[i] = file.name
file.close()
except (ConnectionError, FileNotFoundError):
err = Markdown('### Can not download data')
t_vec, y_ch2 = readwfm(fnames[0])
_, y_ch3 = readwfm(fnames[1])
_, y_ch4 = readwfm(fnames[2])
_, y_ch5 = readwfm(fnames[3])
_, y_ch6 = readwfm(fnames[4])
_, y_ch7 = readwfm(fnames[5])
Check if data has sufficient sample rate
In [6]:
dt = t_vec[1]-t_vec[0]
MSamples = 1e-6/(dt)
downsample = int(np.rint(float(1e-6) / dt))
if MSamples >= 120:
out = Markdown(f'### Data with {MSamples:.1f}MS/s resolution')
low_res_data = False
else:
out = Markdown('### Low resolution for peak recononstruction \n adjust oscilloscope settings')
low_res_data = True
out
Out[6]:
Data with 250.0MS/s resolution¶
In [7]:
# class for storing channels and their parameters and results
@dataclass
class channel:
name : str
time : np.ndarray
y : np.ndarray
noPlot : bool = False
noFit : bool = False
## parameters for peak detection
mean_noise_threshold_mult : float = 1.75
peak_prominence_thr : float = 0.001
# peak rise / fall time
rise_time : float = None
fall_time : float = None
# calibration
calibration : float = None # keV/V
# Flip waveform of PMTs but no SiPMs
flip_waveform : bool = True
## parameters recieved from further analysis
mean_channel_noise : float = 0
# threshold
threshold : float = 0
# arrays for results
simple_peak_loc : np.ndarray = None
lonely_peak_loc : np.ndarray = None
# for debug
lonely_peak_mask : np.ndarray = None
no_corr_peak_height : np.ndarray = None
lonely_peak_height : np.ndarray = None
simple_corr_peak_height : np.ndarray = None
simple_corr_peak_err : np.ndarray = None
simple_corr_peak_filtered_loc : np.ndarray = None
simple_corr_peak_filtered_height : np.ndarray = None
multipass_peak_loc : np.ndarray = None
multipass_corr_peak_height : np.ndarray = None
multipass_corr_peak_err : np.ndarray = None
multipass_corr_peak_filtered_loc : np.ndarray = None
multipass_corr_peak_filtered_height : np.ndarray = None
y_diff : np.ndarray = None
y_gen : np.ndarray = None
y_gen2 : np.ndarray = None
@property
def out_fname(self):
return f'{ch.name}_{shot_no}'
In [8]:
# set channel and their parameters
ch2 = channel('LYSO-3', t_vec, y_ch2, rise_time=6.124e-09, fall_time=5.355e-08, calibration=4884., flip_waveform=False) # @ C:300mV
ch3 = channel('CeBr-a', t_vec, y_ch3, rise_time=1.457e-09, fall_time=2.754e-08, calibration=12578., noFit= True, noPlot = True) # @ HV:600 V
ch4 = channel('CeBr-b', t_vec, y_ch4, rise_time=3.410e-09, fall_time=2.786e-08, calibration=6714.) # @ HV:600 V
ch5 = channel('NaITl', t_vec, y_ch5, rise_time=2.258e-08, fall_time=2.202e-07, calibration=71799., noFit = True) # @ HV:600 V
ch6 = channel('LYSO-1', t_vec, y_ch6, rise_time=5.946e-09, fall_time=5.146e-08, calibration=5007., flip_waveform=False) # @ C:300mV
ch7 = channel('LYSO-2', t_vec, y_ch7, rise_time=5.840e-09, fall_time=4.985e-08, calibration=3713., flip_waveform=False) # @ C:300mV
channels = [ch3,ch4,ch5,ch6,ch7,ch2]
In [9]:
color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axes = plt.subplots(len(channels), sharex =True)
[ax.plot(ch.time[::downsample],ch.y[::downsample],
label = ch.name, color = c) for ax, ch, c in zip(axes, channels, color_cycle)]
[ax.axhline(1., color = 'red', ls = '--')
for ax, ch in zip(axes, channels) if not ch.flip_waveform and np.any(ch.y > 1.)]
fig.legend()
is_clipping = [np.any(ch.y > 1.) for ch in channels if not ch.flip_waveform]
if np.any(is_clipping):
warning = Markdown("""### SiPM's integrated amplifier saturated
HXR spectrum reconstruction from whole discharge may be unreliable
""")
display(warning)
SiPM's integrated amplifier saturated¶
HXR spectrum reconstruction from whole discharge may be unreliable
Remove offset from OSC waveforms and flip¶
In [10]:
def get_noise_level(t,y_vals,noise_time):
y_noise = y_vals[t < noise_time]
y_noise = np.nan_to_num(y_noise, posinf = 0, neginf = 0)
mean_y = np.mean(y_noise)
# mean noise amplitude
abs_noise = np.abs(y_noise - mean_y)**2
mean_noise_ampl = np.sqrt(np.mean(abs_noise))
return mean_y, mean_noise_ampl
def moving_average(x, w):
return np.convolve(x, np.ones(w), 'save') / w
for ch in channels:
noise, ch.mean_channel_noise = get_noise_level(ch.time, ch.y, 0.000)
ch.y = ch.y - noise
ch.threshold = ch.mean_channel_noise*ch.mean_noise_threshold_mult
if ch.flip_waveform:
ch.y *=-1
## smooth-out noise
ch.y = moving_average(ch.y, 5)
Plot OSC waveforms and peak thresholds¶
In [11]:
#%matplotlib widget
fig, axes = plt.subplots(len(channels), sharex =True)
for ax,ch,color in zip(axes,channels, color_cycle):
ax.plot(ch.time[::downsample]*1e3, ch.y[::downsample], label = ch.name, color=color)
ax.axhline(ch.mean_channel_noise,ls= '-.', c = 'k')
ax.axhline(ch.threshold,ls= '--', c = 'k')
ax.set_ylim(-2*ch.threshold, 2*ch.threshold)
ax.grid()
fig.legend();
Find peaks¶
In [12]:
max_HXR_energy = 10.
for ch in channels:
peaks_loc, _ = find_peaks(ch.y, height=ch.threshold, prominence = ch.mean_channel_noise*2)
ch.simple_peak_loc = peaks_loc
ch.no_corr_peak_height = ch.y[peaks_loc]
if ch.calibration is not None and not ch.noPlot and not ch.noFit:
if ch.name == 'NaITl':
continue
max_HXR_energy = np.max((max_HXR_energy, ch.no_corr_peak_height.max()*ch.calibration))
Define useful functions¶
In [13]:
@vectorize(['float64(float64, boolean)'])
def half_exp(x, where):
'''compute exp function only for selected indices (to avoid infinite/large values)'''
if where:
return np.exp(x)
else:
return 0.0
def get_index(oneDArray, val):
diff = np.abs(oneDArray - val)
return np.argmin(diff)
def exp_decay(t, height, center, rise, fall):
'''smooth peak fucntion from https://arxiv.org/pdf/1706.01211.pdf eq:(2)'''
# bool mask where compute exp functions (without it results in +-inf values)
exp_window_split_index = get_index(t, center-100*rise)
exp_window = np.zeros_like(t, np.bool_)
exp_window[exp_window_split_index:].fill(True)
# scaling and shifting of peak fucntion
A = height / ( (fall-rise)/fall * (fall/rise - 1 )**(-rise/fall) )
peak_center_shift = rise*np.log((fall-rise)/rise)
center = center - peak_center_shift
peak = A * 1/(1+half_exp(-(t-center)/rise, exp_window)) * half_exp(-(t-center)/fall, exp_window)
return peak
def gen_peaks(time, HXR_y, peaks_loc, channel_props, test_diff: bool = True):
corrected_peak_heights = np.zeros_like(peaks_loc, float)
HXR_generated = np.zeros_like(HXR_y)
diff = HXR_y.copy()
discarded_peaks = np.zeros_like(peaks_loc, np.bool_)
# operate only on small slice of waveform (it's faster)
dt = np.mean(np.diff(time))
i_slice = max(int(channel_props.fall_time*100 / dt),1000)
for n, peak in enumerate(peaks_loc):
# use differece between generated and original waveform to obtain peak height
peak_height = diff[peak]
peak_center = time[peak]
slice_idx_left = max(peak-i_slice,0)
slice_idx_right = min(peak+i_slice,len(HXR_y)-1)
gen_slice_time = time[slice_idx_left:slice_idx_right]
if peak_height < 0.:
corrected_peak_heights[n] = 0.
continue
# generate peak
generated_peak = exp_decay(gen_slice_time, peak_height, peak_center,
channel_props.rise_time, channel_props.fall_time)
# if adding this peak increases differance btw original and generated signal
# discard this peak as false detection
diff_test = diff[slice_idx_left:slice_idx_right].copy()
diff_test -= generated_peak
if test_diff and np.sum(np.abs(diff_test)) > np.sum(np.abs(diff[slice_idx_left:slice_idx_right])):
corrected_peak_heights[n] = 0.
discarded_peaks[n] = True
else:
# subtract peak from signal and record peak height
diff[slice_idx_left:slice_idx_right] -= generated_peak
HXR_generated[slice_idx_left:slice_idx_right] += generated_peak
corrected_peak_heights[n] = peak_height
return corrected_peak_heights, HXR_generated, diff, discarded_peaks
@njit
def is_lonely_peak(time, HXR_y, peaks, peak_no , fall_time, mean_noise_level):
peak_fall_time = fall_time * np.log(HXR_y[peaks[peak_no]]/mean_noise_level)
if peak_no != 0:
t_before = time[peaks[peak_no]] - peak_fall_time
t_peak_before = time[peaks[peak_no-1]]
if t_before < t_peak_before:
return False
if peak_no <= peaks.size-2:
t_after = time[peaks[peak_no]] + peak_fall_time
t_peak_after = time[peaks[peak_no+1]]
if t_peak_after < t_after:
return False
return True
def identify_lonely_peaks(time,HXR_y, peaks, channel_props):
def is_lonely_peak_reduced(n):
return is_lonely_peak(time,
HXR_y,
peaks,
n,
channel_props.fall_time,
channel_props.mean_channel_noise)
lonely_peaks_idx = [is_lonely_peak_reduced(n) for n in range(peaks.size)]
lonely_peaks_idx = np.array(lonely_peaks_idx, np.bool_)
return lonely_peaks_idx
def find_nearest_idx(array,value):
idx = np.searchsorted(array, value, side="left")
if idx > 0 and (idx == len(array) or math.fabs(value - array[idx-1]) < math.fabs(value - array[idx])):
return idx-1
else:
return idx
def estimate_peak_err(time, HXR_y, diff, peak_loc, channel_props):
fall_time = channel_props.fall_time
mean_noise_level = channel_props.mean_channel_noise
peak_fall_times = fall_time * np.log(HXR_y[peak_loc]/mean_noise_level)
# get maximum of diff in [peak-peak_fall_time ; peak+peak_fall_time]
left_index = (find_nearest_idx(time, time[p]-peak_fall_time) for p,peak_fall_time in zip(peak_loc,peak_fall_times))
right_index = (find_nearest_idx(time, time[p]+peak_fall_time) for p, peak_fall_time in zip(peak_loc,peak_fall_times))
# how to vectorize?
peak_error_estimate = np.array([np.max(np.abs(diff[l:r])) for l,r in zip(left_index,right_index)])
return peak_error_estimate
In [14]:
for ch in channels:
if ch.noFit or ch.rise_time is None or ch.rise_time is None:
continue
ch.simple_corr_peak_height, ch.y_gen, ch.y_diff ,_ = gen_peaks(ch.time, ch.y,
ch.simple_peak_loc,
ch)
ch.simple_corr_peak_err = estimate_peak_err(ch.time, ch.y, ch.y_diff, ch.simple_peak_loc, ch)
ch.lonely_peak_mask = identify_lonely_peaks(ch.time, ch.y,
ch.simple_peak_loc,
ch)
ch.lonely_peak_loc = ch.simple_peak_loc[ch.lonely_peak_mask].copy()
ch.lonely_peak_height = ch.no_corr_peak_height[ch.lonely_peak_mask].copy()
relative_error = np.zeros_like(ch.simple_corr_peak_height)
np.divide(ch.simple_corr_peak_err, ch.simple_corr_peak_height, out=relative_error, where=ch.simple_corr_peak_height != 0)
filtering_mask = relative_error < 0.1
ch.simple_corr_peak_filtered_loc = ch.simple_peak_loc[filtering_mask]
ch.simple_corr_peak_filtered_height = ch.simple_corr_peak_height[filtering_mask]
In [15]:
# plots = list()
# for ch in channels:
# if ch.noFit or ch.rise_time is None or ch.rise_time is None or ch.noPlot:
# continue
# l_gen = '%s gen' % ch.name
# l_dif = '%s diff' % ch.name
# data_to_plot = {ch.name : ch.y, l_gen : ch.y_gen, l_dif : ch.y_diff}
# df = pd.DataFrame(data_to_plot,index = pd.Index(ch.time*1e3, name = 'time'))
# ddf = dd.from_pandas(df, npartitions=os.cpu_count())
# gen_pl = ddf[l_gen].hvplot.line(rasterize = True, cmap=['green'], colorbar=False)
# org_pl = ddf[ch.name].hvplot.line(rasterize = True, cmap=['blue'], colorbar=False)
# diff_pl = ddf[l_dif].hvplot.line(rasterize = True, cmap=['red'], colorbar=False)
# peaks_pl = hv.ErrorBars((ch.time[ch.simple_peak_loc]*1e3,ch.y[ch.simple_peak_loc],ch.simple_corr_peak_height,0),
# vdims=['U_%s' % ch.name, 'yerrneg', 'yerrpos'],label = 'peaks').opts(upper_head=None, lower_head=None)
# plot = org_pl * gen_pl * diff_pl * peaks_pl
# plot.opts(width=950, height = 200).relabel(ch.name)
# plots.append(plot)
# hv.Layout(plots).cols(1)
In [16]:
fig,(ax,ax2,ax3) = plt.subplots(3, sharex = True, sharey=True)
for ch, color in zip(channels, color_cycle):
ax.hist(ch.no_corr_peak_height, bins = 100, label = '%s' % ch.name, histtype = 'step', color = color)
if ch.noFit or ch.rise_time is None or ch.fall_time is None:
continue
ax2.hist(ch.simple_corr_peak_height[ch.lonely_peak_mask], bins = 100, histtype = 'step', color = color)
ax3.hist(ch.simple_corr_peak_filtered_height, bins = 100, histtype = 'step', color = color)
ax.grid()
ax2.grid()
ax3.grid()
ax.semilogy()
ax2.semilogy()
ax3.semilogy()
fig.legend();
Use differance to find less prominent peaks¶
In [17]:
for ch in channels:
if ch.noFit or ch.rise_time is None or ch.rise_time is None:
continue
peaks_loc = ch.simple_peak_loc
for pass_no in range(3):
peaks_locNew, _ = find_peaks(ch.y_diff, height=ch.threshold, prominence = ch.mean_channel_noise*2)
if peaks_locNew.size == 0:
continue
_, _, ch.y_diff,_ = gen_peaks(ch.time,
ch.y_diff,
peaks_locNew,
ch,
test_diff = False)
peaks_loc = np.concatenate((peaks_loc,peaks_locNew))
peaks_loc = np.sort(peaks_loc)
ch.multipass_peak_loc = peaks_loc
Re-generate waveform¶
In [18]:
plots = list()
for ch in channels:
if ch.noFit or ch.rise_time is None or ch.rise_time is None:
continue
peaks_loc = ch.multipass_peak_loc
corrected_peak_heights, ch.y_gen2, ch.y_diff, discarded = gen_peaks(ch.time,
ch.y,
peaks_loc,
ch,
test_diff = True)
### Discard false positives, peaks below threshold and saturated values
saturated = ch.y[peaks_loc]>0.9
mask_discarded = np.logical_or(discarded, corrected_peak_heights<ch.threshold, saturated)
corrected_peak_heights = corrected_peak_heights[~mask_discarded]
peaks_loc = peaks_loc[~mask_discarded]
ch.multipass_peak_loc = peaks_loc
ch.multipass_corr_peak_height = corrected_peak_heights
ch.multipass_corr_peak_err = estimate_peak_err(ch.time, ch.y, ch.y_diff, ch.multipass_peak_loc, ch)
relative_error = np.zeros_like(ch.multipass_corr_peak_height)
np.divide(ch.multipass_corr_peak_err, ch.multipass_corr_peak_height, out=relative_error, where=ch.multipass_corr_peak_height != 0)
filtering_mask = relative_error < 0.1
ch.multipass_corr_peak_filtered_loc = ch.multipass_peak_loc[filtering_mask]
ch.multipass_corr_peak_filtered_height = ch.multipass_corr_peak_height[filtering_mask]
In [19]:
# plots = list()
# for ch in channels:
# if ch.noFit or ch.rise_time is None or ch.rise_time is None or ch.noPlot:
# continue
# l_gen = '%s gen' % ch.name
# l_gen2 = '%s gen2' % ch.name
# l_dif = '%s diff' % ch.name
# data_to_plot = {ch.name : ch.y, l_gen : ch.y_gen, l_dif : ch.y_diff, l_gen2 : ch.y_gen2}
# df = pd.DataFrame(data_to_plot, index = pd.Index(ch.time*1e3, name = 'time'))
# ddf = dd.from_pandas(df, npartitions=os.cpu_count())
# gen_pl = ddf[l_gen].hvplot.line(rasterize = True, cmap=['olive'], colorbar=False)
# gen_pl2 = ddf[l_gen2].hvplot.line(rasterize = True, cmap=['green'], colorbar=False)
# org_pl = ddf[ch.name].hvplot.line(rasterize = True, cmap=['blue'], colorbar=False)
# diff_pl = ddf[l_dif].hvplot.line(rasterize = True, cmap=['red'], colorbar=False)
# peaks_pl = hv.ErrorBars((ch.time[ch.multipass_peak_loc]*1e3,ch.y[ch.multipass_peak_loc],ch.multipass_corr_peak_height,0),
# vdims=['U_%s' % ch.name, 'yerrneg', 'yerrpos']).opts(upper_head=None, lower_head=None)
# plot = org_pl * gen_pl * gen_pl2 * diff_pl * peaks_pl
# plot.opts(width=950, height = 200).relabel(ch.name)
# plots.append(plot)
# plot = hv.Layout(plots).cols(1)
# plot
Show HXR energy spectrum¶
In [20]:
fig, ax = plt.subplots()
for ch,color in zip(channels,color_cycle):
if ch.noPlot or ch.calibration is None:
continue
_, bins, _ = ax.hist(ch.no_corr_peak_height*ch.calibration, bins = 100, label = ch.name,
range = (0,max_HXR_energy), histtype = 'step', ls = '--', color = color)
if ch.noFit or ch.rise_time is None or ch.rise_time is None:
continue
h1, b1 = np.histogram(ch.simple_corr_peak_filtered_height*ch.calibration, bins = len(bins)-1, range = (bins[0], bins[-1]))
h2, b2 = np.histogram(ch.multipass_corr_peak_filtered_height*ch.calibration, bins = len(bins)-1, range = (bins[0], bins[-1]))
h = (h1+h2)/2
h_p = np.maximum(h1,h2)-h
h_m = h-np.minimum(h1,h2)
ax.errorbar(bins[:-1], h, yerr = (h_m, h_p), color = color, ls = '')
fig.legend(frameon =True, fancybox=True,framealpha=1.)
ax.semilogy()
ax.grid()
ax.set_xlabel('$E$ [keV]')
ax.set_ylabel('N [-]')
Out[20]:
Text(0, 0.5, 'N [-]')
In [21]:
fig,(ax,ax2,ax3,ax4) = plt.subplots(4, sharex = True, sharey=True)
for ch, color in zip(channels, color_cycle):
if ch.noFit or ch.rise_time is None or ch.fall_time is None:
continue
ax.hist(ch.simple_corr_peak_height*ch.calibration, bins = 100, label = '%s' % ch.name, histtype = 'step', color = color)
ax2.hist(ch.simple_corr_peak_filtered_height*ch.calibration, bins = 100, histtype = 'step', color = color)
ax3.hist(ch.multipass_corr_peak_height*ch.calibration, bins = 100, histtype = 'step', color = color)
ax4.hist(ch.multipass_corr_peak_filtered_height*ch.calibration, bins = 100, histtype = 'step', color = color)
ax.grid()
ax2.grid()
ax3.grid()
ax4.grid()
ax.semilogy()
ax2.semilogy()
ax3.semilogy()
ax4.semilogy()
fig.legend();
HXR spectrum in different parts of discharge¶
In [22]:
# waveform_plots = list()
# for ch,color in zip(channels,color_cycle):
# if ch.noPlot or ch.calibration is None:
# continue
# w_plot = decimate(hv.Curve((ch.time*1000, ch.y), 't [ms]', 'U [V]', label = ch.name))
# w_plot.opts(width = 600, height = 150)
# waveform_plots.append(w_plot)
# # parameters for sliders
# window = .002 # 1ms
# time_range = np.arange(0, t_vec[-1]*1000,1)
# l_hists_hmap = hv.g (kdims=['Time'])
# s_hists_hmap = hv.HoloMap(kdims=['Time'])
# m_hists_hmap = hv.HoloMap(kdims=['Time'])
# span_hmap = hv.HoloMap(kdims=['Time'])
# for time_sel in time_range:
# time_sel/=1000
# s_hist_plots, m_hist_plots, l_hist_plots = [],[],[]
# left_margin, right_margin = time_sel - window/2, time_sel + window/2
# for ch in channels:
# if ch.noFit or ch.rise_time is None or ch.rise_time is None or ch.calibration is None or ch.noPlot:
# continue
# for peak_loc, peak_height, hist_plots in zip([ch.lonely_peak_loc,ch.simple_peak_loc,ch.multipass_peak_loc],
# [ch.lonely_peak_height,ch.simple_corr_peak_height, ch.multipass_corr_peak_height],
# [l_hist_plots,s_hist_plots,m_hist_plots]):
# peak_time = ch.time[peak_loc]
# mask_peaks = (peak_time > left_margin) & (peak_time < right_margin)
# # histograms for selected time slice
# peak_heights = peak_height[mask_peaks]
# frequencies, edges = np.histogram(peak_heights*ch.calibration, 100, range = (0,max_HXR_energy))
# h_plot = hv.Curve((edges, frequencies), 'E [keV]', 'N [-]', label = ch.name)
# h_plot.opts(logy=True, ylim = (1,None), xlim = (0,None), width = 600, height = 200)
# hist_plots.append(h_plot)
# time_sel*=1000
# l_hists_hmap[time_sel] = hv.Overlay(l_hist_plots)
# s_hists_hmap[time_sel] = hv.Overlay(s_hist_plots)
# m_hists_hmap[time_sel] = hv.Overlay(m_hist_plots)
# span_hmap[time_sel] = hv.VSpan(left_margin*1e3,right_margin*1e3)
# #waveform_plots = hv.Overlay(waveform_plots)
# #hv.Layout((waveform_plots*span_hmap)+ l_hists_hmap+s_hists_hmap+m_hists_hmap).cols(1).opts(title = '')
# hv.Layout(l_hists_hmap+s_hists_hmap+m_hists_hmap).cols(1).opts(title = '')
Peak-count variability in discharge¶
In [23]:
slice_size = 1 # per ms
time_slices = np.arange(0, t_vec[-1]*1e3,slice_size)
LYSO1_pc, LYSO2_pc, LYSO3_pc = np.zeros_like(time_slices), np.zeros_like(time_slices), np.zeros_like(time_slices)
LYSO1_pcM, LYSO2_pcM, LYSO3_pcM = np.zeros_like(time_slices), np.zeros_like(time_slices), np.zeros_like(time_slices)
for time_slice_idx, time_slice in enumerate(time_slices):
for ch, pc, pcM in zip([ch6,ch7,ch2],[LYSO1_pc,LYSO2_pc,LYSO3_pc],[LYSO1_pcM,LYSO2_pcM,LYSO3_pcM]):
peak_times = ch.time[ch.simple_peak_loc]
pc[time_slice_idx] = np.sum(np.logical_and(time_slice < peak_times*1e3, peak_times*1e3 < time_slice + slice_size))
peak_times = ch.time[ch.simple_peak_loc]
pcM[time_slice_idx] = np.sum(np.logical_and(time_slice < peak_times*1e3, peak_times*1e3 < time_slice + slice_size))
fig, ax = plt.subplots(1)
ax.plot(time_slices,LYSO1_pc, label = ch6.name)
ax.plot(time_slices,LYSO2_pc, label = ch7.name)
ax.plot(time_slices,LYSO3_pc, label = ch2.name)
ax2.plot(time_slices,LYSO1_pcM, label = ch6.name)
ax2.plot(time_slices,LYSO2_pcM, label = ch7.name)
ax2.plot(time_slices,LYSO3_pcM, label = ch2.name)
ax.legend()
ax.set_ylabel('#(simple) [-]')
ax2.set_ylabel('#(corr) [-]')
ax2.set_xlabel('t [ms]');
In [24]:
LYSOA_ratio = np.divide(LYSO1_pc, LYSO2_pc, out=np.ones_like(LYSO1_pc), where=LYSO2_pc!=0)
LYSOB_ratio = np.divide(LYSO3_pc, LYSO2_pc, out=np.ones_like(LYSO1_pc), where=LYSO2_pc!=0)
LYSOA_ratioM = np.divide(LYSO1_pcM, LYSO2_pcM, out=np.ones_like(LYSO1_pcM), where=LYSO2_pcM!=0)
LYSOB_ratioM = np.divide(LYSO3_pcM, LYSO2_pcM, out=np.ones_like(LYSO1_pcM), where=LYSO2_pcM!=0)
fig, (ax,ax2) = plt.subplots(2, sharex = True)
ax.set_yscale('log')
ax2.set_yscale('log')
ax.plot(time_slices,LYSOA_ratio, label = f'{ch6.name}/{ch7.name}')
ax.plot(time_slices,LYSOB_ratio, label = f'{ch2.name}/{ch7.name}')
ax2.plot(time_slices,LYSOA_ratioM, label = f'{ch6.name}/{ch7.name}')
ax2.plot(time_slices,LYSOB_ratioM, label = f'{ch2.name}/{ch7.name}')
ax.legend()
ax.set_ylabel('ratio(simple) [-]')
ax2.set_ylabel('ratio(corr) [-]')
ax2.set_xlabel('t [ms]');
In [25]:
LYSO_ratio = np.divide(LYSO1_pc, LYSO3_pc, out=np.ones_like(LYSO1_pc), where=LYSO3_pc!=0)
LYSO_ratioM = np.divide(LYSO1_pcM, LYSO3_pcM, out=np.ones_like(LYSO1_pcM), where=LYSO3_pcM!=0)
fig, (ax,ax2) = plt.subplots(2, sharex = True)
ax.plot(time_slices,LYSOA_ratio, label = f'{ch6.name}/{ch2.name}')
ax2.plot(time_slices,LYSO_ratioM, label = f'{ch6.name}/{ch2.name}')
ax.legend()
ax.set_ylabel('ratio(simple) [-]')
ax2.set_ylabel('ratio(corr) [-]')
ax2.set_xlabel('t [ms]');
ax.set_yscale('log')
ax2.set_yscale('log')
In [26]:
ip_url = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/Results/Ip.csv'
ip = ds.open(ip_url)
ip = np.loadtxt(ip,delimiter=',')
In [27]:
plot_channels = []
for ch,color in zip(channels,color_cycle):
if ch.noPlot or ch.calibration is None or ch.noFit:
continue
plot_channels.append(ch)
fig, axes = plt.subplots(nrows = len(plot_channels), ncols =2, figsize = [1.5*6.4, 1.75*4.8])
prevAx = None
prevrAx = None
for ch,color,ax, rax in zip(plot_channels,color_cycle,axes[:,0],axes[:,1]):
if ch.noPlot or ch.calibration is None:
continue
_, bins, _ = ax.hist(ch.no_corr_peak_height*ch.calibration, bins = 100, label = ch.name,
range = (0,max_HXR_energy), histtype = 'step', ls = '--', color = color)
if ch.noFit or ch.rise_time is None or ch.rise_time is None:
continue
h1, b1 = np.histogram(ch.simple_corr_peak_height*ch.calibration,
bins = len(bins)-1, range = (bins[0], bins[-1]))
h2, b2 = np.histogram(ch.multipass_corr_peak_height*ch.calibration,
bins = len(bins)-1, range = (bins[0], bins[-1]))
h = (h1+h2)/2
h_p = np.maximum(h1,h2)-h
h_m = h-np.minimum(h1,h2)
ax.errorbar(bins[:-1], h, yerr = (h_m, h_p),ds = 'steps', color = color)
h3, b3 = np.histogram(ch.simple_corr_peak_height[ch.lonely_peak_mask]*ch.calibration,
bins = len(bins)-1, range = (bins[0], bins[-1]))
ax.step(bins[:-1], h3, color = color, ls = '--', alpha = .75)
ip_l = rax.plot(ip[:,0],ip[:,1],'k--', label = '$I_p$')
rax.grid(True,axis='x')
rax = rax.twinx()
peak_time = ch.time[ch.simple_peak_loc] * 1e3
peak_height = ch.simple_corr_peak_height*ch.calibration
rax.scatter(peak_time,peak_height, color = color,s = .1 )
rax.grid(True)
#rax.yaxis.tick_right()
if not ch.flip_waveform:
max_reliableE = ch.calibration * .9
if np.any(ch.y[ch.simple_peak_loc]>max_reliableE):
rax.axhline(max_reliableE/2,color = 'r',ls = '-')
ax.semilogy()
ax.grid()
ax.set_ylabel('N [-]')
if prevAx is not None:
ax.sharex(prevAx)
prevAx = ax
rax.yaxis.set_label_position("right")
rax.set_ylabel('$E$ [keV]')
if prevrAx is not None:
#rax.sharex(prevrAx)
rax.sharey(prevrAx)
prevrAx = rax
newHandles = [ip_l[0]]
for ax in axes[:,0]:
handles, labels = ax.get_legend_handles_labels()
newHandles += handles
fig.legend(handles= newHandles,frameon =True, fancybox=True,framealpha=1.,loc = 'lower center')
ax.set_xlabel('$E$ [keV]')
rax.set_xlabel('t [ms]');
rax.set_ylim(-100,max_HXR_energy)
fig.suptitle(f"#{shot_no}");
plt.savefig('icon-fig.png')
