from numpy import *
from matplotlib.pyplot import *

calibration = 662 / 0.641382 #[keV] Cs 137 662 keV  ... X V at 400 V on scintilator
maximum_energy = 1e3 #[keV] maximum to plot in detail
maximum_DAS_voltage = 10 #[V] maximum possible measured voltage
levels = 500 #levels in histogram
omit_peaks = range(25) #number of Be peaks to omit during Cs 137 calibration 

data = loadtxt('HXR_')
peaks = zeros(len(data) - 2) #empty array for peaks
d_data = diff(data) #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] = data[2:][location] #copy only the peaks
peaks *= calibration #calibrate volatge to get energy

counts, bins, patches = hist(peaks, linspace(0.1, maximum_DAS_voltage * calibration, levels)) #plot histogram from 0.1 to max with 1e3 levels

ylabel("Impulse count")
xlabel("Energy [keV]")
savefig("HXR_spectrum_full.png")
xlim(xmax=maximum_energy) #detail plot
savefig("HXR_spectrum_zoom.png")
clf()

counts[omit_peaks]  = 0   #omit the first Be peak

max_indx = argmax( counts )

#saveconst( "max_peak",  )