In [597]:
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 [598]:
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)]))

Constant Intensity

Variables are poorly named and multiply used.


In [599]:
def Poisson(n,nav):
    return (nav**n)*((factorial(n, exact=False))**(-1))*np.exp(-nav)

Rate 1000


In [600]:
ThermR1000 = np.loadtxt('ConstantIntensity_Rate1000.csv',delimiter = ',');

In [601]:
np.max(ThermR1000)


Out[601]:
7.0

In [602]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR1000)+2)

no, binso, patcheso = plt.hist(ThermR1000,normed = 0,align='left',bins = bins,facecolor='green', alpha=0.75);
yerr = np.sqrt(no)/np.sum(no)
plt.xlim(-0.5,np.max(ThermR1000)+1);



In [603]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR1000)+2)

n, bins, patches = plt.hist(ThermR1000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50);

plt.scatter(np.arange(0,np.max(ThermR1000)+1), n)
plt.errorbar(np.arange(0,np.max(ThermR1000)+1),n,yerr,fmt='none')


plt.xlim(-0.5,np.max(ThermR1000)+1);
plt.ylim(0,0.45);


This data was collected in a time window of 1 ms. Need average number of photons counted in 1ms to continue.


In [604]:
n


Out[604]:
array([ 0.273,  0.356,  0.215,  0.104,  0.037,  0.009,  0.004,  0.002])

In [605]:
#n_av = 0*P(0) + 1*P(1) + ....
n_av_CI_R1000 = np.sum(n*np.arange(0,np.max(ThermR1000)+1))
n_av_CI_R1000


Out[605]:
1.329

In [606]:
unc_n_av_CI_R1000 = (1/1000)*np.sqrt(np.sum(ThermR1000))
unc_n_av_CI_R1000


Out[606]:
0.03645545226711637

In [607]:
CHI_CIR1000 = chi2(no,Poisson(np.arange(0,np.max(ThermR1000)+1),n_av_CI_R1000)*1000)
CHI_CIR1000


Out[607]:
0.55076137481786358

In [608]:
CHI_CIR1000/8


Out[608]:
0.068845171852232948

In [609]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR1000)+2)

n, bins, patches = plt.hist(ThermR1000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50,label='Data');

plt.errorbar(np.arange(0,np.max(ThermR1000)+1),n,yerr,fmt='none')

plt.plot(np.linspace(0,np.max(ThermR1000)+1), Poisson(np.linspace(0,np.max(ThermR1000)+1),n_av_CI_R1000), 'r--', linewidth=1,label='Poisson')
plt.legend(loc='best')


plt.xlim(-0.5,np.max(ThermR1000)+1);
plt.ylim(0,0.45);

plt.text(4.7,0.3,'$\chi^2 = %.2f$ ' % (CHI_CIR1000),size=20)
plt.text(4.7,0.25,'$n_{av} = %.2f \pm %.2f$ $ms^{-1}$ ' % (n_av_CI_R1000,unc_n_av_CI_R1000),size=20)


plt.xlabel('Counts per ms',size = 20);
plt.ylabel('Probability Distribution',size = 20);
plt.savefig('CIR1000.png')


Rate 3000


In [610]:
ThermR3000 = np.loadtxt('ConstantIntensity_Rate3000.csv',delimiter = ',');

In [611]:
np.max(ThermR3000)


Out[611]:
12.0

In [612]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR3000)+2)

no, binso, patcheso = plt.hist(ThermR3000,normed = 0,align='left',bins = bins,facecolor='green', alpha=0.75);
yerr = np.sqrt(no)/np.sum(no)
plt.xlim(-0.5,np.max(ThermR3000)+1);



In [613]:
len(no)


Out[613]:
13

In [614]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR3000)+2)

n, bins, patches = plt.hist(ThermR3000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50);

plt.scatter(np.arange(0,np.max(ThermR3000)+1), n)
plt.errorbar(np.arange(0,np.max(ThermR3000)+1),n,yerr,fmt='none')


plt.xlim(-0.5,np.max(ThermR3000)+1);
plt.ylim(0,0.25);



In [615]:
#n_av = 0*P(0) + 1*P(1) + ....
n_av_CI_R3000 = np.sum(n*np.arange(0,np.max(ThermR3000)+1))
n_av_CI_R3000


Out[615]:
3.9649999999999999

In [616]:
unc_n_av_CI_R3000 = (1/1000)*np.sqrt(np.sum(ThermR3000))
unc_n_av_CI_R3000


Out[616]:
0.062968245965724659

In [617]:
CHI_CIR3000 = chi2(no,Poisson(np.arange(0,np.max(ThermR3000)+1),n_av_CI_R3000)*1000)
CHI_CIR3000


Out[617]:
0.6398657508282598

In [618]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR3000)+2)

n, bins, patches = plt.hist(ThermR3000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50,label='Data');
plt.errorbar(np.arange(0,np.max(ThermR3000)+1),n,yerr,fmt='none')
plt.plot(np.linspace(0,np.max(ThermR3000)+1), Poisson(np.linspace(0,np.max(ThermR3000)+1),n_av_CI_R3000), 'r--', linewidth=1,label='Poisson')


plt.legend(loc='best')
plt.xlim(-0.5,np.max(ThermR3000)+1);
plt.ylim(0,0.25);

plt.text(8,0.16,'$\chi^2 = %.2f$ ' % (CHI_CIR3000),size=20)
plt.text(8,0.12,'$n_{av} = %.2f \pm %.2f$ $ms^{-1}$ ' % (n_av_CI_R3000,unc_n_av_CI_R3000),size=20)

plt.xlabel('Counts per ms',size = 20);
plt.ylabel('Probability Distribution',size = 20);
plt.savefig('CIR3000.png')


Rate 10,000


In [619]:
ThermR10000 = np.loadtxt('ConstantIntensity_Rate10000.csv',delimiter = ',');

In [620]:
np.max(ThermR10000)


Out[620]:
24.0

In [621]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR10000)+2)

no, binso, patcheso = plt.hist(ThermR10000,normed = 0,align='left',bins = bins,facecolor='green', alpha=0.75);
yerr = np.sqrt(no)/np.sum(no)
plt.xlim(-0.5,np.max(ThermR10000)+1);



In [622]:
len(no)


Out[622]:
25

In [623]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR10000)+2)

n, bins, patches = plt.hist(ThermR10000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50);

plt.scatter(np.arange(0,np.max(ThermR10000)+1), n)
plt.errorbar(np.arange(0,np.max(ThermR10000)+1),n,yerr,fmt='none')


plt.xlim(-0.5,np.max(ThermR10000)+1);
plt.ylim(0,0.13);



In [624]:
#n_av = 0*P(0) + 1*P(1) + ....
n_av_CI_R10000 = np.sum(n*np.arange(0,np.max(ThermR10000)+1))
n_av_CI_R10000


Out[624]:
10.19

In [625]:
unc_n_av_CI_R10000 = (1/1000)*np.sqrt(np.sum(ThermR10000))
unc_n_av_CI_R10000


Out[625]:
0.10094552986635912

In [626]:
CHI_CIR10000 = chi2(no,Poisson(np.arange(0,np.max(ThermR10000)+1),n_av_CI_R10000)*1000)
CHI_CIR10000


Out[626]:
0.5079198794036397

In [627]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR10000)+2)

n, bins, patches = plt.hist(ThermR10000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50,label='Data');
plt.errorbar(np.arange(0,np.max(ThermR10000)+1),n,yerr,fmt='none')
plt.plot(np.linspace(0,np.max(ThermR10000)+1,200), Poisson(np.linspace(0,np.max(ThermR10000)+1,200),n_av_CI_R10000), 'r--', linewidth=1,label='Poisson')


plt.legend(loc='best')
plt.xlim(-0.5,np.max(ThermR10000)+1);
plt.ylim(0,0.13);

plt.text(15,0.09,'$\chi^2 = %.2f$ ' % (CHI_CIR10000),size=20)
plt.text(15,0.07,'$n_{av} = %.2f \pm %.2f$ $ms^{-1}$ ' % (n_av_CI_R10000,unc_n_av_CI_R10000),size=20)

plt.xlabel('Counts per ms',size = 20);
plt.ylabel('Probability Distribution',size = 20);
plt.savefig('CIR10000.png')


Psuedothermal Source


In [628]:
def BoseEin(n,nav):
    return (nav**n)/((nav + 1)**(n+1))

Rate 1000


In [629]:
ThermR1000 = np.loadtxt('Therm_Rate1000.csv',delimiter = ',');

In [630]:
np.max(ThermR1000)


Out[630]:
12.0

In [631]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR1000)+2)

no, binso, patcheso = plt.hist(ThermR1000,normed = 0,align='left',bins = bins,facecolor='green', alpha=0.75);
yerr = np.sqrt(no)/np.sum(no)
plt.xlim(-0.5,np.max(ThermR1000)+1);



In [632]:
len(no)


Out[632]:
13

In [633]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR1000)+2)

n, bins, patches = plt.hist(ThermR1000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50);

plt.scatter(np.arange(0,np.max(ThermR1000)+1), n)
plt.errorbar(np.arange(0,np.max(ThermR1000)+1),n,yerr,fmt='none')


plt.xlim(-0.5,np.max(ThermR1000)+1);
plt.ylim(0,0.45);



In [634]:
#n_av = 0*P(0) + 1*P(1) + ....
n_av_CI_R1000 = np.sum(n*np.arange(0,np.max(ThermR1000)+1))
n_av_CI_R1000


Out[634]:
1.2669999999999999

In [635]:
unc_n_av_CI_R1000 = (1/1000)*np.sqrt(np.sum(ThermR1000))
unc_n_av_CI_R1000


Out[635]:
0.035594943461115372

In [636]:
CHI_ThermR1000 = chi2(no,BoseEin(np.arange(0,np.max(ThermR1000)+1),n_av_CI_R1000)*1000)
CHI_ThermR1000


Out[636]:
2.1680114946865223

In [637]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR1000)+2)

n, bins, patches = plt.hist(ThermR1000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50,label='Data');
plt.errorbar(np.arange(0,np.max(ThermR1000)+1),n,yerr,fmt='none')
plt.plot(np.linspace(0,np.max(ThermR1000)+1,200), BoseEin(np.linspace(0,np.max(ThermR1000)+1,200),n_av_CI_R1000), 'r--', linewidth=1,label='Bose-Einstein')


plt.legend(loc='best')

plt.xlim(-0.5,np.max(ThermR1000)+1);
plt.ylim(0,0.45);

plt.text(7.2,0.31,'$\chi^2 = %.2f$ ' % (CHI_ThermR1000),size=20)
plt.text(7.2,0.26,'$n_{av} = %.2f \pm %.2f$ $ms^{-1}$ ' % (n_av_CI_R1000,unc_n_av_CI_R1000),size=20)

plt.xlabel('Counts per ms',size = 20);
plt.ylabel('Probability Distribution',size = 20);
plt.savefig('ThermR1000.png')


Rate 3000


In [638]:
ThermR3000 = np.loadtxt('Therm_Rate3000.csv',delimiter = ',');

In [639]:
np.max(ThermR1000)


Out[639]:
12.0

In [640]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR3000)+2)

no, binso, patcheso = plt.hist(ThermR3000,normed = 0,align='left',bins = bins,facecolor='green', alpha=0.75);
yerr = np.sqrt(no)/np.sum(no)
plt.xlim(-0.5,np.max(ThermR3000)+1);



In [641]:
len(no)


Out[641]:
27

In [642]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR3000)+2)

n, bins, patches = plt.hist(ThermR3000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50);

plt.scatter(np.arange(0,np.max(ThermR3000)+1), n)
plt.errorbar(np.arange(0,np.max(ThermR3000)+1),n,yerr,fmt='none')


plt.xlim(-0.5,np.max(ThermR3000)+1);
plt.ylim(0,0.25);



In [643]:
#n_av = 0*P(0) + 1*P(1) + ....
n_av_CI_R3000 = np.sum(n*np.arange(0,np.max(ThermR3000)+1))
n_av_CI_R3000


Out[643]:
3.077

In [644]:
unc_n_av_CI_R3000 = (1/1000)*np.sqrt(np.sum(ThermR3000))
unc_n_av_CI_R3000


Out[644]:
0.055470712993434652

In [645]:
CHI_ThermR3000 = chi2(no,BoseEin(np.arange(0,np.max(ThermR3000)+1),n_av_CI_R3000)*1000)
CHI_ThermR3000


Out[645]:
2.3866423125870884

In [646]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR3000)+2)

n, bins, patches = plt.hist(ThermR3000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50,label='Data');
plt.errorbar(np.arange(0,np.max(ThermR3000)+1),n,yerr,fmt='none')
plt.plot(np.linspace(0,np.max(ThermR3000)+1,200), BoseEin(np.linspace(0,np.max(ThermR3000)+1,200),n_av_CI_R3000), 'r--', linewidth=1,label='Bose-Einstein')


plt.legend(loc='best')

plt.xlim(-0.5,np.max(ThermR3000)+1);
plt.ylim(0,0.25);

plt.text(16,0.17,'$\chi^2 = %.2f$ ' % (CHI_ThermR3000),size=20)
plt.text(16,0.15,'$n_{av} = %.2f \pm %.2f$ $ms^{-1}$ ' % (n_av_CI_R3000,unc_n_av_CI_R3000),size=20)

plt.xlabel('Counts per ms',size = 20);
plt.ylabel('Probability Distribution',size = 20);
plt.savefig('ThermR3000.png')


Rate 10,000


In [647]:
ThermR10000 = np.loadtxt('Therm_Rate10000.csv',delimiter = ',');

In [648]:
np.max(ThermR10000)


Out[648]:
96.0

In [649]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR10000)+2)

no, binso, patcheso = plt.hist(ThermR10000,normed = 0,align='left',bins = bins,facecolor='green', alpha=0.75);
yerr = np.sqrt(no)/np.sum(no)
plt.xlim(-0.5,np.max(ThermR10000)+1);



In [650]:
len(no)


Out[650]:
97

In [651]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR10000)+2)

n, bins, patches = plt.hist(ThermR10000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50);

plt.scatter(np.arange(0,np.max(ThermR10000)+1), n)
plt.errorbar(np.arange(0,np.max(ThermR10000)+1),n,yerr,fmt='none')


plt.xlim(-0.5,np.max(ThermR10000)+1);
plt.ylim(0,0.11);



In [652]:
#n_av = 0*P(0) + 1*P(1) + ....
n_av_CI_R10000 = np.sum(n*np.arange(0,np.max(ThermR10000)+1))
n_av_CI_R10000


Out[652]:
10.186

In [653]:
unc_n_av_CI_R10000 = (1/1000)*np.sqrt(np.sum(ThermR10000))
unc_n_av_CI_R10000


Out[653]:
0.10092571525632107

In [654]:
CHI_ThermR10000 = chi2(no,BoseEin(np.arange(0,np.max(ThermR10000)+1),n_av_CI_R10000)*1000)
CHI_ThermR10000


Out[654]:
1.1986746354337972

In [655]:
plt.figure(figsize=(10,6))
plt.xticks(size = 13);
plt.yticks(size = 13);
bins = np.arange(0,np.max(ThermR10000)+2)

n, bins, patches = plt.hist(ThermR10000,normed = 1,align='left',bins = bins,facecolor='blue', alpha=0.50,label='Data');
plt.errorbar(np.arange(0,np.max(ThermR10000)+1),n,yerr,fmt='none')
plt.plot(np.linspace(0,np.max(ThermR10000)+1,200), BoseEin(np.linspace(0,np.max(ThermR10000)+1,200),n_av_CI_R10000), 'r--', linewidth=1,label='Bose-Einstein')


plt.legend(loc='best')

plt.xlim(-0.5,np.max(ThermR10000)+1);
plt.ylim(0,0.11);

plt.text(57,0.075,'$\chi^2 = %.2f$ ' % (CHI_ThermR10000),size=20)
plt.text(57,0.06,'$n_{av} = %.2f \pm %.2f$ $ms^{-1}$ ' % (n_av_CI_R10000,unc_n_av_CI_R10000),size=20)

plt.xlabel('Counts per ms',size = 20);
plt.ylabel('Probability Distribution',size = 20);
plt.savefig('ThermR10000.png')


Determining the discriminator level


In [656]:
V = np.array([-5,-10,-15,-20,-30,-40,-50,-60,-80,-100,-120]) #Discriminator voltage. 1s Counting
Counts = np.array([658984,549870,459153,346578,138705,28503,3638,889,381,191,130]) #photon counts
BkdCounts = np.array([32.4,23.5,21.9,21.5,18.5,14.2,13,7.5,2.4,0.4,0.2]) #bkd

RealCounts = Counts - BkdCounts #bkd subtracted

S2N = RealCounts/BkdCounts #Signal to noise ratio

In [657]:
fig, ax1 = plt.subplots(figsize=(10,6))

ax1.scatter(V, S2N,c='b',marker="o",label='Signal-to-Noise')
ax1.set_xlabel('Discriminator Voltage (mV)',size=20)
# Make the y-axis label, ticks and tick labels match the line color.
ax1.set_ylabel('Signal-to-Noise', color='b',size=20)
ax1.tick_params('y', colors='b')
plt.legend(loc='lower left')
plt.xlim(-120,0)

ax2 = ax1.twinx()
ax2.scatter(V, RealCounts*1e-5,c='r',marker=".",label='Real Counts')
ax2.set_ylabel('Real Counts $\cdot 10^{5}$', color='r',size=20)
ax2.tick_params('y', colors='r')
plt.xlim(-120,0)
plt.ylim(-1.3,7)


fig.tight_layout()
plt.savefig('Discriminator.png')
plt.legend(loc='lower right')
plt.show()


Pulse Observation


In [658]:
Const_I_1 = pd.read_excel('const_i_1.csv',skiprows=[0]);

Time1 = Const_I_1['second']
Volt1 = Const_I_1['Volt']

plt.figure(figsize=(6,3))
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.plot(Time1*1e6,Volt1*1e3)
plt.ylim(-160,80)
plt.xlim(-5,5);



In [659]:
Const_I_2 = pd.read_excel('const_i_2.csv',skiprows=[0]);

Time2 = Const_I_2['second']
Volt2 = Const_I_2['Volt']

plt.figure(figsize=(6,3))
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.plot(Time2*1e6,Volt2*1e3)
plt.ylim(-110,50)
plt.xlim(-5,5);



In [660]:
Const_I_3 = pd.read_excel('const_i_3.csv',skiprows=[0]);

Time3 = Const_I_3['second']
Volt3 = Const_I_3['Volt']

plt.figure(figsize=(6,3))
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.plot(Time3*1e6,Volt3*1e3)
plt.ylim(-180,40);
plt.xlim(-5,5);



In [661]:
Const_I_4 = pd.read_excel('const_i_4.csv',skiprows=[0]);

Time4 = Const_I_4['second']
Volt4 = Const_I_4['Volt']

plt.figure(figsize=(6,3))
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.plot(Time4*1e6,Volt4*1e3)
plt.ylim(-120,50)
plt.xlim(-5,5);



In [662]:
f, axarr = plt.subplots(2, 2);
axarr[0, 0].plot(Time1*1e6,Volt1*1e3);
#axarr[0, 0].set_xlabel('Time ($\mu$s)');
axarr[0, 0].set_ylabel('Pulse Height (mV)');
axarr[0, 0].set_xlim(-5,5);
axarr[0, 0].set_ylim(-160,80);

axarr[0, 1].plot(Time2*1e6,Volt2*1e3);
#axarr[0, 1].set_xlabel('Time ($\mu$s)');
#axarr[0, 1].set_ylabel('Pulse Height (mV)');
axarr[0, 1].set_xlim(-5,5);
axarr[0, 1].set_ylim(-110,50);

axarr[1, 0].plot(Time3*1e6,Volt3*1e3);
axarr[1, 0].set_xlabel('Time ($\mu$s)');
axarr[1, 0].set_ylabel('Pulse Height (mV)');
axarr[1, 0].set_xlim(-5,5);
axarr[1, 0].set_ylim(-180,40);

axarr[1, 1].plot(Time4*1e6,Volt4*1e3);
axarr[1, 1].set_xlabel('Time ($\mu$s)');
#axarr[1, 1].set_ylabel('Pulse Height (mV)');
axarr[1, 1].set_xlim(-5,5);
axarr[1, 1].set_ylim(-120,50);
# Fine-tune figure; hide x ticks for top plots and y ticks for right plots
fig.tight_layout()
#plt.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False);
#plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=False);
plt.savefig('Pulses.png')


Compare BE with P


In [663]:
def Poisson(n,nav):
    return (nav**n)*((factorial(n, exact=False))**(-1))*np.exp(-nav)

def BoseEin(n,nav):
    return (nav**n)/((nav + 1)**(n+1))

In [668]:
n = np.linspace(0,100)

In [691]:
f, axarr = plt.subplots(1, 2,figsize=(15,4));
n = np.linspace(0,10)


plt.sca(axarr[0])
plt.plot(n,Poisson(n,1),label='Poisson');
plt.plot(n,BoseEin(n,1),label='Bose-Einstein');
#axarr[0, 0].set_xlabel('Time ($\mu$s)');
plt.ylabel('Probability Distribution');
plt.text(4,0.22,'$n_{av} = 1$',size=20)
plt.xlabel('Counts per unit time');
plt.xlim(0,6)
#axarr[0, 0].set_xlim(-5,5);
#axarr[0, 0].set_ylim(-160,80);
plt.legend()

plt.sca(axarr[1])
n = np.linspace(0,20)
plt.plot(n,Poisson(n,4),label='Poisson');
plt.plot(n,BoseEin(n,4),label='Bose-Einstein');
plt.xlim(0,12)
plt.xlabel('Counts per unit time');
plt.text(8,0.11,'$n_{av} = 4$',size=20)
#axarr[0, 1].set_xlabel('Time ($\mu$s)');
#axarr[0, 1].set_ylabel('Pulse Height (mV)');
#axarr[0, 1].set_xlim(-5,5);
#axarr[0, 1].set_ylim(-110,50);
plt.legend()


plt.savefig('BEPCompare.png')



In [ ]: