In [89]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit # import the curve fitting function
from scipy.special import factorial
%matplotlib inline
In [90]:
#Lets determine scale from channel calibration
x = np.array([94, 186, 279, 370, 464]) #Channels
x = x-20 #Measurement issues
t = np.array([2, 4, 6, 8, 10]) #Delay in us
#Fitting
def myfun(x,a,b):
ans = a*x + b
return ans
p0 = [0.216,0.4] #guess
xspace = np.linspace(0,500,1000) #FIX?
plsq, pcov = curve_fit(myfun, x, t, p0) # curve fit returns p and covariance matrix
# these give the parameters and the uncertainties
a= plsq[0]
ea = np.sqrt(pcov[0,0])
b = plsq[1] #microseconds
eb = np.sqrt(pcov[1,1]) #microseconds
yfit = myfun(xspace,plsq[0],plsq[1]) # use fit results for a, b, c
plt.figure(figsize=(10,6));
#plt.scatter(Channels,Counts);
plt.scatter(x,t,color='red',label='Data');
plt.plot(xspace,yfit,label='Fit')
plt.legend(loc='upper left')
plt.xlabel('Channel Number',fontsize=20);
plt.ylabel('Delay Time [$\mu s$]',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.xlim(0,500)
plt.ylim(0,12)
plt.text(290,4,'$a = %.4f \pm %.4f$ $\mu s / N$' % (a,ea),size=20)
plt.text(290,3,'$b = %.2f \pm %.2f$ $\mu s$ ' % (b,eb),size=20)
plt.savefig('Calibration.png')
In [91]:
Day1 = pd.read_excel('DecayData.xlsx',sheetname = None, header = None)
Channels = np.arange(0,512) #Channel Number
Counts = np.array([entry for entry in Day1['Sheet1'][0]]) #Counts
#We want to group Counts
N = 5
GroupedChannels = np.arange(0,512,N)
GroupedCounts = np.array([np.mean(Counts[i:i+N]) for i in GroupedChannels])
#(1/N)*SQRT(Sum of Squares of Individual Uncertainties)
GroupedUncertainties = np.array([(1/N)*np.sqrt(np.sum(Counts[i:i+N])) for i in GroupedChannels])
GroupedTimes = GroupedChannels*a + b #Need to determine Scale above [Scale] = microseconds/channel
GroupedCountsALL = GroupedCounts[2:]
GroupedUncertaintiesALL = GroupedUncertainties[2:]
GroupedTimesALL = GroupedTimes[2:]
#Lop off first two entries
I = 15
F = -1
GroupedCounts = GroupedCounts[I:F]
GroupedUncertainties = GroupedUncertainties[I:F]
GroupedTimes = GroupedTimes[I:F]
#Fitting
def myfun(t,No,tau,c):
ans = No*np.exp(-t/tau) + c
return ans
p0 = [20,2,0.1] #guess
xspace = np.linspace(0,12,1000)
plsq, pcov = curve_fit(myfun, GroupedTimes, GroupedCounts, p0) # curve fit returns p and covariance matrix
# these give the parameters and the uncertainties
No= plsq[0]
eNo = np.sqrt(pcov[0,0])
tau = plsq[1] #microseconds
etau = np.sqrt(pcov[1,1]) #microseconds
c = plsq[2]
ec = np.sqrt(pcov[2,2])
yfit = myfun(xspace,plsq[0],plsq[1],plsq[2]) # use fit results for a, b, c
plt.figure(figsize=(10,6));
#plt.scatter(Channels,Counts);
plt.plot(xspace,yfit)
plt.scatter(GroupedTimesALL,GroupedCountsALL,color='red');
plt.errorbar(GroupedTimesALL,GroupedCountsALL, yerr= GroupedUncertaintiesALL,fmt='None')
plt.xlabel('Delay Time [$\mu s$]',fontsize=20);
plt.ylabel('Counts',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.xlim(0,12)
plt.ylim(0,17)
plt.text(7,13,'$tau = %.2f \pm %.2f$ $\mu s$ ' % (tau,etau),size=20)
plt.savefig('Day1.png')
In [92]:
def chi2(ymeasured,ytheory):
# N = number of bins
# ymeasured = Array of the measured number of counts in each bin
# ytheory = Array of the predicted number of counts in each bin
N = len(ymeasured)
sigma = np.sqrt(ymeasured)
for i in np.arange(0,N):
if sigma[i] == 0:
sigma[i] = 1
return (1/N)*np.sum(np.array([((ymeasured[i]-ytheory[i])/sigma[i])**2 for i in np.arange(0,N)]))
In [96]:
Counts.shape
Out[96]:
In [95]:
Day2 = pd.read_excel('ArrivalData.xlsx',sheetname = None, header = None)
#Read in
Channels = np.arange(0,512) #Channel Number
Counts = np.array([entry for entry in Day2['Sheet1'][0]]) #Counts
#Initiate Hist
plt.figure(figsize=(10,6));
n, bins, patches = plt.hist(Counts,bins=11,align='left',label='Data')
#Uncertainties
nerr = np.sqrt(n)
plt.errorbar(np.arange(1,12), n, yerr = nerr ,fmt='None')
#Poisson Distribution
def poisson(k,mu):
return np.exp(-mu) * mu**k / factorial(k,0)
mu = np.mean(Counts)
y = poisson(bins, mu)
l = plt.plot(bins, y*512, 'r--', linewidth=1,label='Poisson')
#Chi Squared
no = np.sum(bins[0:11]*n)/np.sum(n)
X2 = chi2(n,y*512)
plt.text(7.5,90,'$\chi^2 = %.2f$ ' % (X2),size=20)
plt.text(7.5,82,'$n_{av} = %.2f \pm %.2f$ $s^{-1}$ ' % (no,np.sqrt(no)),size=20)
plt.xlabel('Events in One Second',fontsize=20);
plt.ylabel('Occurence Number',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.xlim(0,12);
plt.legend(loc='upper left')
plt.savefig('Day2.png')
In [ ]: