I'm currently in the process of computing a differential cross-section for the scattering of a 420 MeV electron by an O-16 nucleus (with a Wood-Saxon charge distribution). The problem is that the resulting graphs don't look anything like they should. I first calculated the form factor for the wood-saxon charge distribution (which is spherically symmetric) using the quad function of scipy.integrate then calculated the Rutherford cross-section, the recoil factor which are quite straightforward. The differential cross-section is then given by
$$\frac{d\sigma}{d\Omega} = \left(\frac{d\sigma}{d\Omega}\right)_\text{Ruth}\left(\frac{1}{1 + \frac{q^2}{2ME}}\right)|F(q)|^2$$
With $M$ the mass of the nucleus, $E$ the energy of the incoming electron, $q$ the transfered momentum and
$$F(q) = \frac{4\pi}{q}\int_0^\infty dr \frac{r\sin(qr/\hbar)}{1+e^{(r-R)/a}}$$
I used natural units for the different quantities. I suspect the problem to come from the form factor but how exactly I don't know. I tried to improve the integration accuracy as much as I could but it didn't help. Also, somehow, having a smaller step for the angle array seems to make it worse. Here is the Python code I made :
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from scipy.special import expit
from scipy.signal import savgol_filter
def wood_saxon_model(R, a, r):
density = expit((R - r) / a)
return density
def form_factor(q):
R = 1.07 * 16*(1/3)
a = 0.5
#hbar = 6.582119e-22
hbar = 1
integral_result, _ = integrate.quad(lambda r: r np.sin(q * r / hbar) * wood_saxon_model(R, a, r), 0, np.inf, limit = 1000)
return 4 * np.pi * hbar * integral_result / q
def rutherford_cs(Z1, Z2, q):
alpha = 7.29735256e-3
#hbarc = 1.97327e2
hbarc = 1
cs = (Z1 * Z2 * alpha * hbarc)2 / (q4 + 1e-10)
return cs
def recoil(E, M, q):
return 1 / (1 + q*2 / (2 M * E))
def cross_section(Z1, Z2, E, M, q):
return rutherford_cs(Z1, Z2, q) * recoil(E, M, q) * form_factor(q)**2
theta = np.linspace(np.pi / 18, 6 * np.pi / 18, 1000)
Z1 = 1
Z2 = 8
E = 420
M = 1.489916863e4
q = 2 * E * np.sin(theta / 2)
wood_saxon_cross_section = [cross_section(Z1, Z2, E, M, x) for x in q]
squared_form_factor = [form_factor(x)**2 for x in q]
#smoothed_cross_section = savgol_filter(wood_saxon_cross_section, window_length=11, polyorder=3)
plt.plot(np.degrees(theta), squared_form_factor)
plt.show()
plt.plot(np.degrees(theta), wood_saxon_cross_section)
plt.yscale('log')
plt.xlabel('Scattering Angle (degrees)')
plt.ylabel('Differential Cross Section')
plt.title('Differential Cross Section')
plt.show()
I'm clueless as to what might be the trouble here. Here is a comparison of what I get vs what kind of thing I should get:
What I have vs what I should have (image source)
