plt.figure(figsize=(10,6));
plt.scatter(Peaks,Energy);
plt.xlim(0,240)
plt.ylim(0,1000)
plt.xlabel('x (mm)');
plt.ylabel('y (mm)');
plt.plot(xlots,yfit);
plt.legend(['data','Fit'],loc='lower right');
plt.text(5,900,'a = %.3f +/- %.3f keV' % (plsq[0], np.sqrt(pcov[0,0])),size=17)
plt.text(5,800,'b = %.3f +/- %.3f keV MCA$^{-1}$' % (plsq[1], np.sqrt(pcov[1,1])),size=17)
plt.text(5,700,'c = (%.4f +/- %.4f)$\cdot 10^{-3}$ keV MCA$^{-2}$' % (plsq[2]*1e3, np.sqrt(pcov[2,2])*1e3),size=17)
plt.xlabel('MCA Number',fontsize=20);
plt.ylabel('Energy (keV)',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.savefig('LinearMCAFit')

def deltaE(N,dN):
    daN = np.sqrt((da/a)**2 + (dN/N)**2)*(a*N)

    dbN2 = np.sqrt((db/b)**2 + 4*(dN/N)**2)*(b*N**2)

    dcN3 = np.sqrt((dc/c**2) + 9*(dN/N)**2)*(c*N**3)

    dEMeas = np.sqrt(daN**2 + dbN2**2 + dcN3**2)*1e-3 #Convert to KeV
    
    return dEMeas

N = np.array([102.40]) #Channel number of 60 degree scattered photopeak
dN = np.array([9.00])
theta = np.array([60])*np.pi/180

EMeas = myfun(N,a,b,c)
EMeas

dEMeas = deltaE(N,dN)
dEMeas

Eo = 661.7 #Initial keV energy of gamma rays (before scattering) 
mc2 = 511 #electron mass in keV

def ECompton(Eo,mc2,theta):
    return Eo/(1+(Eo/mc2)*(1-np.cos(theta)))

EComp = ECompton(Eo,mc2,theta)
EComp

"""
Explanation: 1
End of explanation
"""


thetas = np.linspace(-np.pi,np.pi,50);

plt.figure(figsize=(10,6));
plt.plot(thetas,ECompton(Eo,mc2,thetas),label='Compton');
plt.errorbar(theta,EMeas,dEMeas);
plt.scatter(theta,EMeas,dEMeas,label='Measured');
plt.legend();

plt.xlabel('Scattering Angle [Radians]',fontsize=20);
plt.ylabel('Final Energy (keV)',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.xlim(-np.pi,np.pi);
#plt.savefig('Sample')

def Thomson(theta):
    ro = 2.82*1e-15
    return (1/2)*(ro**2)*(1+np.cos(theta)**2) #set b = 1

def KleinNishina(theta):
    ro = 2.82*1e-15
    gamma = Eo/mc2
    
    return (1/2)*(ro**2)*(1+np.cos(theta)**2)*((1+gamma*(1-np.cos(theta)))**(-2))*(1+((gamma*(1-np.cos(theta)))**2)/((1+np.cos(theta)**2)*(1+gamma*(1-np.cos(theta)))))

thetas = np.linspace(-np.pi,np.pi,50);

plt.figure(figsize=(10,6));
plt.plot(thetas,Thomson(thetas),label='Thomson');
plt.plot(thetas,KleinNishina(thetas),label='Klein-Nishina');
plt.legend();

plt.xlabel('Scattering Angle [Radians]',fontsize=20);
plt.ylabel('Differential Cross section',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.xlim(-np.pi,np.pi);
#plt.savefig('Sample')

"""
Explanation: 2

Test how the energy of scattered atoms varies with scattering angle
Use (1) to determine keV mass of electron
Which of Thomson and Klein-Nishina differential cross section is a better description?
End of explanation
"""


