# coding: utf-8
# In[1]:
"""
Created on Thu Feb 08 19:50:25 2018
@author: Ondrej
"""
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 14 18:52:48 2015
@author: Ondřej Ficker
Creating HXR histogram from calibration data, data stored in binnary file,
reasonable sampling frequency is 1 MHz.
"""
import sys
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
path = 'C:\\Users\\PC\\Documents\\Data_Compass\\ch0_'
num_file = np.arange(16577,16580,step=1)
#U=200
#CsE=662
#mscint=0.2
#def movingaverage(interval, window_size):
# window= np.ones(int(window_size))/float(window_size)
# return np.convolve(interval, window, 'same')
#fig, ax = plt.subplots()
#fig2,ax2=plt.subplots()
#Y=np.fromfile(path + str(num_file[3]) + file,dtype='>f8', count=-1, sep='')
#k=662/0.05
##setting time axis
#print(np.shape(Y))
#
#X=np.linspace(0,len(Y),len(Y))
#
#ax2.plot(X,movingaverage(Y,5))
#plt.show()
#def fit_func(x,a,b,c,d):
# return a*np.exp(-(x-b)**2/c**2) + d
#
#
#def peakdet(v, delta, x=None):
# """
# Converted from MATLAB script at http://billauer.co.il/peakdet.html
# Returns two arrays
# function [maxtab, mintab]=peakdet(v, delta, x)
# %PEAKDET Detect peaks in a vector
# % [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
# % maxima and minima ("peaks") in the vector V.
# % MAXTAB and MINTAB consists of two columns. Column 1
# % contains indices in V, and column 2 the found values.
# %
# % With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
# % in MAXTAB and MINTAB are replaced with the corresponding
# % X-values.
# %
# % A point is considered a maximum peak if it has the maximal
# % value, and was preceded (to the left) by a value lower by
# % DELTA.
# % Eli Billauer, 3.4.05 (Explicitly not copyrighted).
# % This function is released to the public domain; Any use is allowed.
# """
# maxtab = []
# mintab = []
# if x is None:
# x = np.arange(len(v))
# v = np.asarray(v)
# if len(v) != len(x):
# sys.exit('Input vectors v and x must have same length')
# if not np.isscalar(delta):
# sys.exit('Input argument delta must be a scalar')
# if delta <= 0:
# sys.exit('Input argument delta must be positive')
# mn, mx = np.Inf, -np.Inf
# mnpos, mxpos = np.NaN,np.NaN
# lookformax = True
# for i in np.arange(len(v)):
# this = v[i]
# if this > mx:
# mx = this
# mxpos = x[i]
# if this < mn:
# mn = this
# mnpos = x[i]
# if lookformax:
# if this < mx-delta:
# maxtab.append((mxpos, mx))
# mn = this
# mnpos = x[i]
# lookformax = False
# else:
# if this > mn+delta:
# mintab.append((mnpos, mn))
# mx = this
# mxpos = x[i]
# lookformax = True
#
# return np.array(maxtab), np.array(mintab)
#
## Finding peaks
##maxtab, mintab=peakdet(Y,0.001,X)
##ax.set_ylabel('Counts [-]')
##ax.set_xlabel(r'E [keV]')
##print(len(maxtab))
##ax.hist(maxtab[:,1],np.arange(0.0,0.08,0.001)) #histogram in voltage
##ax.hist(maxtab[:,1],np.linspace(250,1600,80)) #histogram in Energy
#
## saving and showing the figure
##plt.savefig('Hist_to_GOLwiki_cal.png')
#res = []
#for item in num_file:
# print(item)
# Y=np.fromfile(path + file + str(item) + '.bin',dtype='>f8', count=-1, sep='')
# X = np.linspace(0,len(Y),len(Y))
# maxtab, mintab = peakdet(Y,1,X)
# print(len(maxtab))
# res.append(maxtab)
#all_res = np.concatenate(res)
# #histogram in voltage
#plt.hist(all_res[:,1],np.arange(0.0,3,0.01))
#plt.xlabel(r'$U \ [\mathrm{V}]$')
#plt.ylabel(r'$N \ [\mathrm{-}]$')
#
#hist,edges = np.histogram(all_res[:,1],bins=150,range=(0.6,3))
#plt_edges = (edges[:-1]+edges[1:])/2
#plt.plot(plt_edges,hist,'.')
#
#x_min = 2.1
#x_max = 3
##hist_fit = hist[(np.where(plt_edges >= x_min)) and (np.where(plt_edges <= x_max))]
#hist_fit = hist[np.where(plt_edges >= x_min)]
#edges_fit = plt_edges[np.where(plt_edges >= x_min)]
#
#popt,pcov = curve_fit(fit_func,edges_fit,hist_fit)
#
#print('Estimated position of the peak: ' + str(popt[1]))
#print('Error of esmtimation: ' + str(np.sqrt(np.diag(pcov))[1]))
#
#x = np.linspace(0.01,3,100)
#plt.figure()
#plt.plot(plt_edges,hist,'.',label='Data')
#plt.plot(edges_fit,hist_fit,'.',label='Fitted data')
#plt.plot(x,fit_func(x,*popt),label='Fit')
#plt.legend()
#plt.xlabel(r'$U \ [\mathrm{V}]$')
#plt.ylabel(r'$N \ [\mathrm{-}]$')
#
#
#plt.figure()
#plt.hist(all_res[:,1],np.arange(0.0,3,0.02))
#plt.plot(x,fit_func(x,*popt),label='Fit')
#plt.xlabel(r'$U \ [\mathrm{V}]$')
#plt.ylabel(r'$N \ [\mathrm{-}]$')
#plt.xlim([1,3])
#plt.title(r'$U_{Source} = 1100 \ \mathrm{V}$')
#plt.legend()
#
#def fit_func2(x,a,b,c):
# return a*np.exp(b*(x-c))
#
#x0 = 2.81065e+6
#x1 = 2.8107e+6
#index1 = (X > x0)
#index2 = (X < x1)
#index = index1 * index2
#popt,pcov = curve_fit(fit_func2,X[index],Y[index])
#plt.plot(X[index],Y[index])
data = []
for i,item in enumerate(num_file):
Y=np.fromfile(path + str(item),dtype='>f8', count=-1, sep='')
data.append([item,Y])
time = np.linspace(952,1112,len(Y))
plt.figure()
for i in range(len(data)):
plt.plot(time,data[i][1],label='Shot: ' + str(data[i][0]))
plt.legend()
plt.xlabel(r'$t \ [\mathrm{ms}]$')
plt.xlabel(r'$HXR \ [\mathrm{a.u.}] $')
# In[ ]:
