import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from scipy.signal import convolve
from scipy.optimize import curve_fit
import h5py
import tcspcfit
# matplotlib.use('TkAgg')
# file = h5py.File('IRF SPCM 680nm.h5', 'r')
[docs]
file = h5py.File('../IRF 680nm.h5', 'r')
# file = h5py.File('IRF_AQR_EMN.h5', 'r')
# file = h5py.File('IRF.h5', 'r')
# file = h5py.File('IRF NIM.h5', 'r')
# file = h5py.File('beads_both_cards.h5', 'r')
# file = h5py.File('IRF FAST 680nm.h5', 'r')
[docs]
datadict = file['Particle 1']
microtimes = datadict['Micro Times (s)']
[docs]
microtimes = microtimes - np.min(microtimes)
[docs]
differences = np.diff(np.sort(microtimes[:]))
[docs]
channelwidth = np.unique(differences)[1]
t = np.arange(0, 25, channelwidth)
irf, t = np.histogram(microtimes, bins=t)
[docs]
halfmaxind = np.where(irf > 0.5*irf.max())
[docs]
fwhm = t[halfmaxind[0][-1]] - t[halfmaxind[0][0]]
# print(halfmaxind[-1])
# print(fwhm)
# plt.plot(t, irf)
# plt.show()
# file1 = h5py.File('IRF SPCM 680nm.h5', 'r')
# datadict1 = file1['Particle 1']
# microtimes1 = datadict['Micro Times (ns)']
# microtimes1 = microtimes1 - np.min(microtimes1)
# t1 = np.arange(0, 25, channelwidth)
# irf1, t1 = np.histogram(microtimes1, bins=t1)
#
[docs]
centre = t[np.argmax(irf)]
gauss = 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(- (t - centre) ** 2 / (2 * sigma ** 2))
[docs]
gauss = gauss * np.trapz(irf, t)
# # irf = irf + 10 # bg of 1
#
# plt.figure()
# plt.plot(t, gauss)
# plt.plot(t, irf)
# plt.show()
#
# fastamp = 0.1
#
#
[docs]
def fitfunc1(t, tau1):
model = np.exp(-t / tau1)
irs = tcspcfit.colorshift(irf, 20)
convd = convolve(irs, model)
# convd = model
convd = convd * (100 / convd.max()) # peak of 100 is about what you get at low excitation power
return convd
[docs]
def fitfunc(t, amp1, tau1, tau2):
model = amp1 * np.exp(-t / tau1) + (1 - amp1) * np.exp(-t / tau2)
irs = tcspcfit.colorshift(irf, 0.1)
convd = convolve(irs, model)
convd = convd * (100 / convd.max()) # peak of 100 is about what you get at low excitation power
return convd
# convd = fitfunc1(t,1.3)
convd = fitfunc(t,0.5, 0.08, 3.4)
# tau = [[1.3, 0.1, 2.9, 0]]
# convd = fitfunc(t,0.5, 0.08, 3.4)
[docs]
tau = [[0.08, 0.05, 0.1, 0], [3.4, 3.2, 3.6, 0]]
# amp = [[0.5, 0.4, 0.6, 0]]
[docs]
amp = [[0.5, 0.4, 0.6, 0], [0.5, 0, 1, 0]]
[docs]
shift = [0, -300, 300, 0]
[docs]
convdnoise = np.random.poisson(convd)
# convdnoise[convdnoise <= 0] = 0
# convdsum = convdnoise.sum()
# convdnoise = convdnoise / convdsum
# convdnoise = convd
[docs]
bg_est = tcspcfit.FluoFit.estimate_bg(convdnoise)
# print(bg_est * convdsum)
print(bg_est)
# fit = tcspcfit.OneExp(irf, convdnoise, t, channelwidth, tau=tau, shift=shift, bg=[5, 0, 10, 0], method='ml',
# boundaries=[500, 1300, 'Manual', False])
# fit = tcspcfit.OneExp(irf, convdnoise, t, channelwidth, tau=tau, shift=shift, method='ls',
# boundaries=[500, 1300, 'Manual', False])
[docs]
fit = tcspcfit.TwoExp(irf, convdnoise, t, channelwidth, tau=tau, shift=shift, bg=[5, 0, 10, 0], amp=amp, method='ml',
boundaries=[500, 1300, 'Manual', False])
print('# photons: ', convdnoise.sum())
print('# fit photos: ', fit.convd.sum())
print('Tau: ', fit.tau)
print('Amp: ', fit.amp)
print('IRFbg: ', fit.irfbg)
print('bg: ', fit.bg)
print('shift: ', fit.shift)
print('stds: ', fit.stds)
[docs]
resid = (fit.convd - fit.measured) #/ np.sqrt(np.abs(fit.measured))
[docs]
acf = np.correlate(resid, resid, 'full')
[docs]
chisq = np.sum(resid ** 2)
print(chisq)
# plt.figure()
# plt.plot(resid, '.')
# plt.show()
plt.figure()
plt.plot(convdnoise[500:1300])
plt.plot(convd[500:1300])
plt.plot(fit.convd)
# plt.plot(irf)
plt.show()
plt.figure()
plt.plot(fit.measured)
plt.plot(fit.convd)
# plt.plot(fit.irf)
plt.show()
#
# # param, pcov = curve_fit(fitfunc, t, convdnoise, p0=[0.01, 0.08, 3.4])
# # print(param)
# # plt.plot(fitfunc(t, param[0], param[1], param[2]))
# # plt.plot(convdnoise)
# # plt.show()
#
# taus = np.array([])
# for i in range(100):
# print(i)
# convd = fitfunc1(t,1.3)
# tau = [[1.3, 0.1, 2.9, 0]]
# shift = [20, -300, 300, 0]
# convd = convd + 2
# convdsum = convd.sum()
# convdnoise = np.random.poisson(convd)
#
# bg_est = tcspcfit.FluoFit.estimate_bg(convdnoise)
#
# fit = tcspcfit.OneExp(irf, convdnoise, t, channelwidth, tau=tau, shift=shift, bg=[5, 0, 10, 0], method='ml',
# boundaries=[500, 1300, 'Manual', False])
# taus = np.append(taus, fit.tau)
#
# print('tau std = ', np.std(taus))
# plt.figure()
# plt.hist(taus)
# plt.show()