import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
import os
"from scipy.signal import correlate"
from multipletau import autocorrelate, correlate
Load of the data from the two detectors
mat1 = np.loadtxt(r"./Data/FIFO/allineamento_diodo2_effettivo_m1.asc")
mat2 = np.loadtxt(r"./Data/FIFO/allineamento_diodo2_effettivo_m2.asc")
# Tempi in s
dati1 = 1/1.05e8 * mat1[:,0] + 2.44e-12 * mat1[:,1]
dati2 = 1/1.05e8 * mat2[:,0] + 2.44e-12 * mat2[:,1]
dati2 += 35e-6
The acquisition was 3 minutes long. We define 36 windows of 5 seconds each. Inside each window we compute the correlation function in intervals of 1us long. The normalized correlation function is defined as \begin{equation} c\left\{av\right\}\left[k\right] = \dfrac{\sum_n a\left[n+k\right]\,v\left[n\right]}{\left<a\right>\left<v\right>} \end{equation}
%%time
corrVect = []
for i in range(36):
# Time windows of 5s each
timeEdges = np.arange(5*i, 5*(i+1), 1e-6)
binc = timeEdges[:-1] + (timeEdges[1] - timeEdges[0])/2
# Select and histogram in each 5s window the counts in 1us subwindows
cond1 = (dati1 > 5*i) & (dati1 < 5*(i+1) )
cond2 = (dati2 > 5*i) & (dati2 < 5*(i+1) )
h1, _ = np.histogram(dati1[cond1], bins = timeEdges)
h2, _ = np.histogram(dati2[cond2], bins = timeEdges)
m1 = dati1[cond1].mean()
m2 = dati2[cond2].mean()
myCorr = correlate(h1, h2, m=16, deltat=1e-6, normalize = True)
corrVect.append(myCorr.copy())
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:412: RuntimeWarning: divide by zero encountered in double_scalars
if np.abs(traceavg1) / np.median(np.abs(v)) < ZERO_CUTOFF:
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:414: RuntimeWarning: divide by zero encountered in double_scalars
if np.abs(traceavg2) / np.median(np.abs(a)) < ZERO_CUTOFF:
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:435: InvalidMWarning: Input dtype is not float; casting to np.float_!
warnings.warn("Input dtype is not float; casting to np.float_!",
CPU times: total: 55.6 s Wall time: 56.5 s
Then i can mediate over the windows and obtain the mean correlation function (and the error, estimated as the error of the mean (std/sqrt(N)))
corrVect = np.array(corrVect)
corrVect.shape
(36, 161, 2)
xVect = corrVect[0,:,0]
meanCorr = np.mean(corrVect[:,:,1], axis = 0)
meanCorr.shape
(161,)
stdCorr = np.std(corrVect[:,:,1], axis = 0)
stdCorr.shape
(161,)
Let's plot the 36 cross correlation functions
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
for j, i in enumerate(corrVect):
ax.plot(i[:,0], i[:,1], "*", label = f"{j}-th windows")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
The last window's histograms
fig, ax = plt.subplots(2,1)
ax[0].plot(binc, h1, ds = "steps-mid", c = "tab:green")
ax[1].plot(binc, h1, ds = "steps-mid", c = "tab:green")
for a in ax:
a.grid()
plt.show()
%matplotlib inline
def myExp(x, N, td, s, T, tt):
return (1/N) * (1/(1+ x/td)) * (1+ s*x/td)**(-.5) * (1 - T + T*np.exp(-x/tt)) / (1-T)
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1]
p0 = [6.93565896e-01, 50e-6, 6.99182585e+00, 7.56745114e-02,1]
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVect, meanCorr,
sigma = stdCorr, absolute_sigma = True,
maxfev = 1000000, p0=p0, )
#bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) ) )
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g")
ax.plot(xVect, myExp(xVect, *popt), ls = ":", c = "k")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"T:\t{popt[3]}")
print(f"tt:\t{popt[4]}")
N: 0.6351440998561837 td: 0.00026046794056623114 s: 12.545714859148859 T: 0.014269501664092736 tt: 33890.960772910854
# Removing triplet
%matplotlib inline
def myExp(x, N, td, ):
return (1/N) * (1/(1+ x/td)) * (1+ x/(td*25) )**(-.5)
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1]
p0 = [6.93565896e-01, 50e-6, ]
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVect, meanCorr,
sigma = stdCorr, absolute_sigma = True,
maxfev = 1000000, p0=p0, )
#bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) ) )
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g")
ax.plot(xVect, myExp(xVect, *popt), ls = ":", c = "k")
ax.plot(xVect, myExp(xVect, *p0), ls = "--", c = "r")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
N: 0.643489866308663 td: 4.628606638669725e-05
%matplotlib inline
def myExp(x, N, td,s, c ):
return (1/N) * (1/(1+ x/td)) * (1+ x/(td*(s**2)) )**(-.5) +c
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1]
p0 = [6.93565896e-01, 50e-6, 5, 0 ]
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVect, meanCorr,
sigma = stdCorr, absolute_sigma = True,
maxfev = 1000000, p0=p0, )
#bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) ) )
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g", label = "Data")
ax.plot(xVect, myExp(xVect, *popt), ls = ":", c = "k", label = "Fit curve")
ax.plot(xVect, myExp(xVect, *p0), ls = "--", c = "r", label = "Initialized curve")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
ax.legend(fontsize = 14)
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"offset:\t{popt[3]}")
N: 0.6423324213598477 td: 4.595612681603429e-05 s: 5.158288647258937 offset: -0.00011697209462548527
# Conti
N = popt[0]
tauD = popt[1]
s = popt[2]
# Tabulato
D = 4.25e-6 * 1e-4
# Calcoli
r0 = np.sqrt(4 * D * tauD)
print(f"r0 = {r0*1e9:.2f} nm")
# Formula: s = (r0/z0)^2
#s = z0/r0
z0 = r0 *s
print(f"z0 = {z0*1e9:.2f} nm")
vol = np.pi ** (3/2) * z0 * (r0**2)
print(f"volume confocale = {vol*1e3*1e15:.2f} fL")
# Formula: 1/N = 1/veff*C ==> N = Veff*c
Na = 6.02e23
c = N/(vol*1e3*Na)
print(f"Concentrazione: {c*1e9:.2f} nM")
print(f"N: {N}\tvol: {vol}\tNa: {Na}\t")
r0 = 279.51 nm z0 = 1441.79 nm volume confocale = 0.63 fL Concentrazione: 1.70 nM N: 0.6423324213598477 vol: 6.272186460305031e-19 Na: 6.02e+23
# eq completa
%matplotlib inline
def myExp(x, N, td, s, T, tt, c):
return (1/N) * (1/(1+ x/td)) * (1+ x/(td*s**2))**(-.5) * (1 - T + T*np.exp(-x/tt)) / (1-T) +c
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1]
p0 = [6.93565896e-01, 50e-6, 5, 7.56745114e-02,1,0]
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVect, meanCorr,
sigma = stdCorr, absolute_sigma = True,
maxfev = 1000000, p0=p0, )
#bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) ) )
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g", label = "Data")
ax.plot(xVect, myExp(xVect, *popt), ls = ":", c = "k", label = "Fit curve")
ax.plot(xVect, myExp(xVect, *p0), ls = "--", c = "r", label = "Initialized curve")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"T:\t{popt[3]}")
print(f"tt:\t{popt[4]}")
print(f"offset:\t{popt[5]}")
N: 0.7437281331389194 td: 4.5956222795341446e-05 s: 5.158244566878252 T: 0.1363343431472159 tt: 22594.22323630056 offset: -0.0001169715453250161
%matplotlib inline
def myExp(x, N, td, s, T, tt):
return (1/N) * (1/(1+ x/td)) * (1+ s*x/td)**(-.5) * (1 - T + T*np.exp(-x/tt)) / (1-T)
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1]
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVect, meanCorr,
sigma = stdCorr, absolute_sigma = True,
maxfev = 1000000, p0=p0,
bounds = ( (0,0,0,0,1), (1,1,10,1,10) )
)
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g")
ax.plot(xVect, myExp(xVect, *popt), ls = ":", c = "k")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"T:\t{popt[3]}")
print(f"tt:\t{popt[4]}")
N: 0.6407807098849005 td: 0.0002383402046573087 s: 9.999999999999998 T: 3.9769530776476936e-10 tt: 9.715292412118874
%matplotlib inline
def myExp(x, N, td, s, c):
return (1/N) * (1/(1+ x/td)) * (1+ s*x/td)**(-.5) + c
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1]
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVect, meanCorr,
sigma = stdCorr, absolute_sigma = True,
maxfev = 1000000, p0=p0, )
#bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) ) )
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g")
ax.plot(xVect, myExp(xVect, *popt), ls = ":", c = "k")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"offset:\t{popt[3]}")
N: 0.6257904232296319 td: 0.00026162285156846416 s: 12.646909640816657 offset: -9.711286421261734e-05
%matplotlib inline
def myExp(x, N, td, s, T, tt, c):
return (1/N) * (1/(1+ x/td)) * (1+ s*x/td)**(-.5) * (1 - T + T*np.exp(-x/tt)) / (1-T) + c
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1, 0.01]
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVect, meanCorr,
sigma = stdCorr, absolute_sigma = True,
maxfev = 1000000, p0=p0, bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) )
)
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g")
ax.plot(xVect, myExp(xVect, *popt), ls = ":", c = "k")
#ax.legend()
ax.set_title("Normalized cross correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized cross correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"T:\t{popt[3]}")
print(f"tt:\t{popt[4]}")
print(f"Offset:\t{popt[5]}")
N: 0.640780730982292 td: 0.00023834019466660952 s: 9.999999999999934 T: 4.3716459579970013e-08 tt: 9.954124064766036 Offset: 7.424947399001792e-13
On the same data (for now only detector 1) i perform auto correlation, defined as the cross correlation where $a$ and $v$ are the same vector
%%time
corrVectAuto = []
for i in range(36):
# Time windows of 5s each
timeEdges = np.arange(5*i, 5*(i+1), 1e-6)
binc = timeEdges[:-1] + (timeEdges[1] - timeEdges[0])/2
# Select and histogram in each 5s window the counts in 1us subwindows
cond1 = (dati1 > 5*i) & (dati1 < 5*(i+1) )
h1, _ = np.histogram(dati1[cond1], bins = timeEdges)
m1 = dati1[cond1].mean()
myCorrAuto = correlate(h1, h1, m=16, deltat=1e-6, normalize = True)
corrVectAuto.append(myCorrAuto.copy())
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:412: RuntimeWarning: divide by zero encountered in double_scalars
if np.abs(traceavg1) / np.median(np.abs(v)) < ZERO_CUTOFF:
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:414: RuntimeWarning: divide by zero encountered in double_scalars
if np.abs(traceavg2) / np.median(np.abs(a)) < ZERO_CUTOFF:
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:435: InvalidMWarning: Input dtype is not float; casting to np.float_!
warnings.warn("Input dtype is not float; casting to np.float_!",
CPU times: total: 48.4 s Wall time: 50.8 s
corrVectAuto = np.array(corrVectAuto)
xVectAuto = corrVectAuto[0,:,0]
meanCorrAuto = np.mean(corrVectAuto[:,:,1], axis = 0)
stdCorrAuto = np.std(corrVectAuto[:,:,1], axis = 0)
Plot and fit with eq. 4
%matplotlib inline
def myExp(x, N, td, s, T, tt, c):
return (1/N) * (1/(1+ x/td)) * (1+ s*x/td)**(-.5) * (1 - T + T*np.exp(-x/tt)) / (1-T) + c
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1, 0.01]
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVectAuto, meanCorrAuto,
sigma = stdCorrAuto, absolute_sigma = True,
maxfev = 1000000, p0=p0, bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) )
)
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVectAuto, meanCorrAuto, "*g")
ax.plot(xVectAuto, myExp(xVectAuto, *popt), ls = ":", c = "k")
#ax.legend()
ax.set_title("Normalized auto correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized auto correlation function", fontsize = 14)
ax.set_ylim((0,2))
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"T:\t{popt[3]}")
print(f"tt:\t{popt[4]}")
print(f"Offset:\t{popt[5]}")
N: 0.5996429513294332 td: 0.0001858665281973763 s: 9.999999999999996 T: 5.208863118139618e-08 tt: 9.99999999992172 Offset: 1.8396988049125845e-12
and with eq. 3
%matplotlib inline
def myExp(x, N, td, s, c):
return (1/N) * (1/(1+ x/td)) * (1+ s*x/td)**(-.5) + c
# p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 7.56745114e-02,1]
p0 = [6.93565896e-01, 1.95028590e-04, 6.99182585e+00, 1]
popt, pcov = curve_fit(myExp, xVectAuto, meanCorrAuto,
sigma = stdCorrAuto, absolute_sigma = True,
maxfev = 1000000, p0=p0, )
#bounds = ( (0,0,0,0,1,0), (1,1,10,1,10,1) ) )
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g")
ax.plot(xVectAuto, myExp(xVectAuto, *popt), ls = ":", c = "k")
ax.set_ylim((0,2))
#ax.legend()
ax.set_title("Normalized auto correlation function in 5s windows", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized auto correlation function", fontsize = 14)
ax.set_xscale("log")
ax.grid()
plt.show()
print(f"N:\t{popt[0]}")
print(f"td:\t{popt[1]}")
print(f"s:\t{popt[2]}")
print(f"offset:\t{popt[3]}")
N: 0.0024421161267835943 td: 0.0006821845650827018 s: 6514692.188705808 offset: -0.0002442783709074924
mat1 = np.loadtxt(r"./Data/FIFO/allineamento_diodo2_effettivo_m1.asc")
mat2 = np.loadtxt(r"./Data/FIFO/allineamento_diodo2_effettivo_m2.asc")
# Tempi in s
dati1 = 1/1.05e8 * mat1[:,0] + 2.44e-12 * mat1[:,1]
dati2 = 1/1.05e8 * mat2[:,0] + 2.44e-12 * mat2[:,1]
%%time
corrVect = []
for i in range(36):
# Time windows of 5s each
timeEdges = np.arange(5*i, 5*(i+1), 1e-6)
binc = timeEdges[:-1] + (timeEdges[1] - timeEdges[0])/2
# Select and histogram in each 5s window the counts in 1us subwindows
cond1 = (dati1 > 5*i) & (dati1 < 5*(i+1) )
cond2 = (dati2 > 5*i) & (dati2 < 5*(i+1) )
h1, _ = np.histogram(dati1[cond1], bins = timeEdges)
h2, _ = np.histogram(dati2[cond2], bins = timeEdges)
m1 = dati1[cond1].mean()
m2 = dati2[cond2].mean()
myCorr = correlate(h1, h2, m=16, deltat=1e-6, normalize = True)
corrVect.append(myCorr.copy())
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:412: RuntimeWarning: divide by zero encountered in double_scalars
if np.abs(traceavg1) / np.median(np.abs(v)) < ZERO_CUTOFF:
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:414: RuntimeWarning: divide by zero encountered in double_scalars
if np.abs(traceavg2) / np.median(np.abs(a)) < ZERO_CUTOFF:
C:\Users\steca\anaconda3\lib\site-packages\multipletau\core.py:435: InvalidMWarning: Input dtype is not float; casting to np.float_!
warnings.warn("Input dtype is not float; casting to np.float_!",
CPU times: total: 56.3 s Wall time: 59 s
corrVect = np.array(corrVect)
corrVect.shape
(36, 161, 2)
xVect = corrVect[0,:,0]
meanCorr = np.mean(corrVect[:,:,1], axis = 0)
meanCorr.shape
(161,)
stdCorr = np.std(corrVect[:,:,1], axis = 0)
stdCorr.shape
(161,)
fig, ax = plt.subplots()
fig.set_size_inches(12,5)
ax.plot(xVect, meanCorr, "*g", label = "NCCF")
#ax.legend()
ax.set_title("NCCF if not shifting the second vector", fontsize = 16)
ax.set_xlabel("Time [s]", fontsize = 14)
ax.set_ylabel("Normalized auto correlation function", fontsize = 14)
ax.axvline(x=35e-6, ls = "--", c = "k", label = "35 us")
ax.set_ylim((0,2))
ax.set_xscale("log")
ax.grid()
ax.legend(fontsize = 14)
plt.show()
