Revision 0ab53102c02cd5c97fc58c71489748c2546cd368 (click the page title to view the current version)
# -*- 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()
