#!/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) 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 def prepare_data(): data = Shot() if data.exist('NIturbo6251'): tvec, signal_0 = data['NIturbo6251'] else: raise IOError('Missing NIturbo6251') start = data['plasma_start'] end = data['plasma_end'] plasma = data['plasma'] 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)])) 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 #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 end = data['plasma_end']*1e3 else: start= 0; end = 40 out = [ [get_data('HXR', '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', '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()