#!/usr/bin/python2
# -*- coding: utf-8 -*-
import matplotlib
matplotlib.rcParams['backend'] = 'Agg'
matplotlib.rc('font', size='10')
matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!!
#import pygolem_lite
from scipy.signal import fftconvolve
from numpy import *
#from pygolem_lite.config import *
from pygolem_lite.modules import *
from matplotlib.pyplot import *
import time
#from shutil import copy, move
import os, sys
from pygolem_lite import Shot
#calibration = 662 / 0.891186 #[keV] Cs 137 662 keV ... X V at 400 V on scintilator
#calibration = 662 / 0.841186 #[keV] Cs 137 662 keV ... X V at 400 V on scintilator (podle mě to blo blbě => o 5% zvýšeno)
#calibration = 662 / 1.3# [keV] Cs 137 662 keV ... novy Vrbuv scintilator
calibration = 662 / 1.15# [keV] Cs 137 662 keV ... X V at 800 V, different amplifier
maximum_energy = 1e3 #[keV] maximum to plot in detail
maximum_DAS_voltage = 10 #[V] maximum possible measured voltage
bins = 500 #bins in histogram
calibration_surface = 5.57 #keV/ms
#min_reliable_intensity = 1000
min_reliable_intensity = 0.5 * calibration # above 0.5 V
def prepare_data():
data = Shot()
if data.exist('nistandard'):
print data['tektronix3014']
tvec, signal_0 = data['tektronix3014'] # 3. channel
signal_0 = signal_0[:,-1]
#tvec, signal_0 = data['ScintProbeNaITl']
#tvec, signals = data['HXR']
if -amin(signal_0) > amax(signal_0):
signal_0 *= -1
else:
raise IOError('Missing nistandard')
start = data['plasma_start']
end = data['plasma_end']
plasma = data['plasma']
maximum_DAS_voltage = amax(signal_0)
signal_0[signal_0 >= maximum_DAS_voltage] = nan
t = time.time()
win = 100
lam = 0.1
signal, chi2 = GapsFilling(signal_0,win,lam) # recontruct cutoffs in the signal_0
print 'filling time', time.time()-t
save_adv('HXR', tvec, signal)
save_adv('HXR_orig', tvec, signal_0)
if plasma:
saveconst('HXR_mean', mean(signal[(tvec > start) & (tvec < end)]))
else:
saveconst('HXR_mean', nan )
sigma = len(tvec)/200 # !!! smoothing parameter !!!!
ker = linspace(-4,4, sigma)
ker = exp( - ker**2)
ker /= sum(ker)
signal[signal < 0 ] = 0
signal_smooth = fftconvolve(signal, ker, mode="same")
signal_smooth*= calibration_surface #calculate detected HXR power (not the total HXR power)
save_adv('HXR_smooth', tvec, signal_smooth)
omit_peaks = range(25) #low intensity noise
peaks = zeros(len(signal) - 2) #empty array for peaks
d_data = diff(signal) #first derivative
d2_data = diff(d_data) #second derivative
location = (d_data[1:] * d2_data < 0) & (d2_data < 0) #location of data where there are peaks
peaks[location] = signal[2:][location] #copy only the peaks
peaks *= calibration #calibrate volatage to get energy
#import IPython
#IPython.embed()
ind_peaks = where(peaks > min_reliable_intensity)[0]
if len(ind_peaks)>1:
first_peak = tvec[ind_peaks[0] ]
saveconst('first_peak', first_peak)
last_peak = tvec[ind_peaks[-1] ]
saveconst('last_peak', last_peak)
#plot histogram from 0.1 to max with 1e3 bins
peaks = peaks[peaks > 0] # remove negative peaks
savetxt('peaks', peaks)
hist, edges = histogram(peaks, bins = bins, range = [0.1, maximum_DAS_voltage * calibration])
#counts, bins, patches = hist(peaks, )
#ylabel("Impulse count")
#xlabel("Energy [keV]")
#savefig("HXR_spectrum_full.png")
#xlim(xmax=maximum_energy) #detail plot
#savefig("HXR_spectrum_zoom.png")
#clf()
hist[omit_peaks] = 0 #omit the first Be peak
max_indx = argmax( hist )
max_peak = edges[1:][max_indx] / calibration
saveconst( "max_peak", max_peak )
print "done"
def graphs(file_type):
data = Shot()
if data['plasma']:
start = data['plasma_start']*1e3-1
end = data['plasma_end']*1e3+1
#last_peak = loadconst('last_peak')*1e3
#end = max(end, last_peak)
else:
start= 0; end = 40
#print "end", end, last_peak, data['plasma_end']*1e3
out = [
[get_data('hxr_signal', 'HXR recalculated', "HXR [a.u.]", xlim=[start, end]),
get_data('HXR_orig', 'HXR', "HXR [a.u.]", xlim=[start, end])],
get_data('HXR_smooth', 'HXR smooth', "HXR power [keV/ms]" , xlim=[start, end], ylim=[0,None]),
get_data('loop_voltage', "Loop voltage" , 'U[V]', xlim=[start, end], ylim=[0,None])
]
print 'multiplot'
multiplot(out, 'Hard-X ray signal' , 'graph', (9,6), 100, 'vertical', file_type)
# icona
os.system('convert -resize 150x120\! graph.png icon.png')
clf()
print '1. hist'
try:
peaks = loadtxt('peaks')
bins = sqrt(len(peaks)) # empirical rule
except:
print "No Xrays were deteced in the signal"
return
#my "intelligent" historogram -- it will make a smaller bins if there are higher points
#density and larger bins if there are only a few points to keep statistical error on the constant value
#fig = figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k')
peaks.sort()
x = hstack((peaks[::int(bins)],peaks[-1]))
y = bins/(diff(x)+1)
#TODO vylapšená verze, otestovat včerně renormalizace
#x = (val[bins/2:]+val[:-(bins)/2])/2
#y = bins/((val[bins/2:]-val[:-(bins)/2]))/2
#y/= sum(diff(x)*y[:-1])
#y[0] = 0
fig = plt.figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k')
semilogx(x[:-1], y, drawstyle = 'steps-post')
axvline(17, color='r', label=r"Molybden K$\alpha$")
axvline(20, color='g', label=r"Molybden K$\beta$")
axvline(662, color='g', label=r"Cesium 662 KeV")
legend(loc = 'best')
ylabel("Impulse count/dE [1/keV]")
xlabel("Energy [keV]")
xlim(10,maximum_DAS_voltage*calibration)
savefig("HXR_spectrum_full.png", bbox_inches='tight', dpi= 100)
clf()
#print '2. hist'
#fig = figure(num=None, figsize=(10, 3), dpi=80, facecolor='w', edgecolor='k')
#bins = sqrt(sum(peaks < maximum_energy)) # empirical rule
#hist_range = linspace(0.1, maximum_energy)
#hist(peaks, hist_range)
#axvline(17, color='r', label=r"Molybden K$\alpha$")
#axvline(20, color='g', label=r"Molybden K$\beta")
#axvline(662, color='g', label=r"Cesium 662 KeV")
#legend(loc = 'best')
#ylabel("Impulse count")
#xlabel("Energy [keV]")
#savefig("HXR_spectrum_zoom.png", bbox_inches='tight', dpi= 100)
#clf()
#print '3. hist'
#scipy.signal.parzen
def postanalysis(file_type):
data = Shot()
if data['plasma']:
start = data['plasma_start']*1e3
end = data['plasma_end']*1e3
else:
start= 0; end = 40
out = [
[get_data('HXR_signal', 'HXR recalculated', "HXR [a.u.]", xlim=[start, end]),
get_data('HXR_orig', 'HXR', "HXR [a.u.]", xlim=[start, end])],
get_data('HXR_smooth', 'HXR smooth', "HXR power [keV/ms]" , xlim=[start, end], ylim=[0,None]),
get_data('loop_voltage', "Loop voltage" , 'U[V]', xlim=[start, end], ylim=[0,None])
]
print 'multiplot'
multiplot(out, 'Hard-X ray signal' , 'graph', (9,6), 100, 'vertical', file_type)
def main():
if sys.argv[1] == "acquisition":
prepare_data()
if sys.argv[1] == "plots":
graphs('png')
saveconst('status', 0)
if sys.argv[1] == "postanalysis":
print "plotting svgz "
postanalysis('svgz')
if __name__ == "__main__":
main()
