# -*- 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

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()
Y=np.fromfile('bin_golem',dtype='>f8', count=-1, sep='')
k=662/1.05
#setting time axis
X=np.linspace(0,len(Y)*(1/1000.0),len(Y))

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(k*Y,0.1,X)
ax.set_ylabel('Counts [-]')
ax.set_xlabel(r'E [keV]')
#ax.hist(maxtab[:,1],np.arange(0.1,3.0,0.01))  #histogram in voltage
ax.hist(maxtab[:,1],np.arange(250,1600,5.8))   #histogram in Energy

# saving and showing the figure
plt.savefig('Hist_to_GOLwiki_cal.png') 
plt.show()