Source code :: main

#!/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'):
	tvec, signals = data['nistandard']
	signal_0 = -signals[:,5] # 6. channel
    else:
	raise IOError('Missing nistandard')
    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)]))
    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

    first_peak = tvec[where(peaks > min_reliable_intensity)[0][0]]
    saveconst('first_peak', first_peak)
    last_peak = tvec[where(peaks > min_reliable_intensity)[0][-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()