# 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[ ]: