/
lifetime-deconvolve
executable file
·109 lines (92 loc) · 3.18 KB
/
lifetime-deconvolve
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#!/usr/bin/env python3
import numpy as np
from numpy import sum, newaxis, sqrt
import scipy.signal, scipy.optimize
from argparse import ArgumentParser
parser = ArgumentParser()
parser.add_argument('corr', metavar='FILE',
help="correlation function")
parser.add_argument('--irf', '-i', metavar='FILE',
help='instrument response function')
parser.add_argument('--components', '-c', type=int, default=1,
help='number of fit components')
parser.add_argument('--rel-resid', '-r', action='store_true',
help='plot relative residuals')
args = parser.parse_args()
jiffy = 4e-12
rep_rate = 80e6
window = 3
print()
print(' %s with IRF from %s' % (args.corr, args.irf))
n = int(window / rep_rate / jiffy)
corr = np.genfromtxt(args.corr)[:n,1]
irf = np.genfromtxt(args.irf)[:n,1]
corr = np.roll(corr, -100)
irf = np.roll(irf, -100)
hi,deconv = scipy.signal.deconvolve(corr, irf)
per = int(1 / rep_rate / jiffy)
times = (np.arange(corr.shape[0]) % per)
# Filter out badness
good_times = times[deconv > 1e-12]
good_corr = deconv[deconv > 1e-12]
# Find beginning of decay
peak_time = np.argmax(deconv) % per
print('Peak time = %f ns' % (peak_time * jiffy / 1e-9))
fit_times = good_times[good_times > peak_time]
fit_corr = good_corr[good_times > peak_time]
# Do fit
def fit_func(t, *components):
offset = components[-1]
a = 0
for i in range(len(components)/2):
tau, amp = components[2*i:2*i+2]
a += amp * np.exp(-t/tau)
return a + offset
p0 = []
for i in range(args.components):
p0.extend([10**i * 1e-9/jiffy, 1e4 / 10**i])
p0.append(0.1)
print(p0)
popt, pcov = scipy.optimize.curve_fit(fit_func, fit_times, fit_corr, p0,
sigma=sqrt(fit_corr))
# Print fit parameters
print('offset = %g' % popt[-1])
for i in range(args.components):
print('component %i' % i)
print(' amplitude = %f' % popt[2*i+1])
print(' tau = %f ns' % (popt[2*i] * jiffy / 1e-9))
print('covariance = ')
print(pcov)
# Plot
from matplotlib import pyplot as pl
pl.suptitle('%s with IRF from %s' % (args.corr, args.irf))
pl.subplot(311)
pl.plot(jiffy * times / 1e-9, corr, '+', label='signal')
pl.plot(jiffy * times / 1e-9, irf, '+', label='IRF')
pl.yscale('log')
pl.legend()
pl.ylabel('counts')
pl.subplot(312)
pl.yscale('log')
pl.plot(jiffy*good_times / 1e-9, good_corr, '+', label='deconvolved')
ts = np.linspace(peak_time, times[-1])
label = ', '.join(['\\tau_%d = %1.3g' % (i, popt[2*i] * jiffy / 1e-9)
for i in range(args.components)])
pl.plot(jiffy*ts / 1e-9, fit_func(ts, *popt), 'r',
label='fit ($%s$ ns)' % label)
pl.legend()
pl.ylabel('counts')
pl.subplot(313)
ts = good_times[good_times > peak_time]
cs = good_corr[good_times > peak_time]
if args.rel_resid:
pl.plot(jiffy * ts / 1e-9, (fit_func(ts, *popt) - cs) / cs, 'b+')
pl.ylim(-0.5, 0.5)
pl.ylabel('rel. residual')
else:
pl.plot(jiffy * ts / 1e-9, fit_func(ts, *popt) - cs, 'b+')
pl.ylabel('residual')
pl.axhline(0, color='k')
#pl.xlim(jiffy*peak_time / 1e-9, jiffy*times[-1] / 1e-9)
pl.xlabel('Lag (ns)')
pl.savefig('%s.png' % args.corr)