#
# -------------------------------
#
# (c) Aurelie Vancraeyenest 2019
# --------------------------------
#
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
[docs]def plotFitHist(data, figname='fig.pdf', bins=161, rmin=2000, rmax=6000):
"""
Plot an histogramm and fit with Gaussian and Lorentzian functions
Parameters
----------
data : np.array
Input data
figname : str
Output filename of the file to save the figure
bins : int
Number of equal-width bins in the given range
rmin : int
max (most right-end bin edge) of the histogramm
rmax : int
min (most left-end bin edge) of the histogramm
See Also
--------
plotHist, plotHists, saveHists
"""
histo, bin_edges = np.histogram(data, bins=bins, range=[rmin, rmax])
bin_centers = 0.5*(bin_edges[1:] + bin_edges[:-1])
# Finding the guess parameters for Gaussian fitting
x0_guess = bin_edges[np.argmax(histo)]
sig_guess = np.std(histo)*4
amp_guess = np.max(histo)
paramL, covL = curve_fit(lorentz, bin_centers, histo)
paramG, covG = curve_fit(gaussian, bin_centers, histo,
p0=[x0_guess, sig_guess, amp_guess])
x_plot = np.linspace(rmin, rmax, 20000)
plt.plot(bin_centers, histo, 'r-')
plt.plot(x_plot, lorentz(x_plot, *paramL), ':m')
plt.plot(x_plot, gaussian(x_plot, *paramG), '--b')
plt.title("60Co resolution")
plt.xlabel("Time")
plt.ylabel("Counts")
plt.legend(('data', 'Fit lorentz', 'Fit Gaussian'),
loc='upper right', shadow=True)
plt.text(rmin, 0.8*np.max(histo),
'Fit parameters:\n FWHM Lor = {:.4}\n FWHM Gaus = {:.5}'
.format(paramL[1], paramG[1]*2.35482))
plt.savefig(figname)
plt.show()
[docs]def plotHist(data, figname='fig.pdf', bins=161, rmin=2000, rmax=6000,
logY=False):
"""
Plot an histogramm and save the output spectrum as a figure
Parameters
----------
data : np.array_like
Input data
figname : str
Output filename of the file to save the figure
bins : int
Number of equal-width bins in the given range
rmin : int
max (most right-end bin edge) of the histogramm
rmax : int
min (most left-end bin edge) of the histogramm
logY : bool
Allow to plot the histogram in a logarithmique scale
See Also
--------
plotHist, plotHists, saveHists
"""
histo, bin_edges = np.histogram(data, bins=bins, range=[rmin, rmax])
bin_centers = 0.5*(bin_edges[1:] + bin_edges[:-1])
if logY is True:
plt.semilogy(bin_centers, histo, 'b-')
else:
plt.plot(bin_centers, histo, 'r-')
plt.title("Time spectrum")
plt.xlabel("Time")
plt.ylabel("Counts")
plt.legend(("data"), loc='upper right', shadow=True)
plt.savefig(figname)
plt.show()
[docs]def plotHists(data, figname='fig.pdf', bins=161, rmin=2000, rmax=6000,
logY=False, style='-'):
"""
Plot multiple histograms and save the output spectrum as a figure
Parameters
----------
data : np.array_like
Input data as multiple numpy arrays
figname : str
Output filename of the file to save the figure
bins : int
Number of equal-width bins in the given range
rmin : int
max (most right-end bin edge) of the histogramm
rmax : int
min (most left-end bin edge) of the histogramm
logY : bool
Allow to plot the histogram in a logarithmique scale
style : str
Plotting style, see matplotlib for more info
See Also
--------
plotHist, plotHists, saveHists
"""
for data in data:
histo, bin_edges = np.histogram(data, bins=bins, range=[rmin, rmax])
bin_centers = 0.5*(bin_edges[1:] + bin_edges[:-1])
if logY is True:
plt.semilogy(bin_centers, histo, '-')
else:
plt.plot(bin_centers, histo, style)
plt.title("Time spectrum")
plt.xlabel("Time")
plt.ylabel("Counts")
plt.legend(("data01", "data02"), loc='upper right', shadow=True)
plt.savefig(figname)
plt.show()
[docs]def saveHists(data, filebase='histo', bins=161, rmin=2000, rmax=6000,
useChnEnding=False):
"""
Sort each data-set in data into histogram and save as text files
Parameters
----------
data : list of np.array_like
Input data asa list of several numpy arrays
filebase : str
Output filename base that will serve to form the output
filename of each text file generated.
Output file format is ASCII
bins : int
Number of equal-width bins in the given range
rmin : int
max (most right-end bin edge) of the histogramm
rmax : int
min (most left-end bin edge) of the histogramm
useChnEnding : bool
If Ture, channel ending type are inserted in the filename
before the file extension. In False, an integer reprenseting
the position in the input dataset is inserted instead
See Also
--------
plotHist, plotHists, saveHists
"""
ii = 0
chnEndings = ('01', '02', '12')
for data in data:
histo, bin_edges = np.histogram(np.array(data), bins=bins,
range=[rmin, rmax])
if useChnEnding is True:
filename = filebase + '_' + chnEndings[ii] + '.hst'
else:
filename = filebase + '_' + str(ii) + '.txt'
np.savetxt(filename, histo, fmt='%i',
header="bins: {}\nrange: [{}:{}]".format(bins, rmin, rmax))
ii += 1
[docs]def lorentz(x, x0, sig, amp):
""" Formula for a Lorentzian
"""
return (amp/(2*np.pi)*sig/((x-x0)**2+sig**2/4))
[docs]def gaussian(x, x0, sig, amp):
""" Formula for a Gaussian
"""
return (amp/(sig*np.sqrt(2*np.pi))*np.exp(-((x-x0)**2)/(2*sig**2)))