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 625.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 = 3.5
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)
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')
