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 / (q**4 + 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.plot(np.degrees(theta), wood_saxon_cross_section)
plt.xlabel('Scattering Angle (degrees)')
plt.ylabel('Differential Cross Section')
plt.title('Differential Cross Section')

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 What I should have

What I have vs what I should have (image source)

  • 2
    $\begingroup$ It seems like you are off by a factor of $10^{30}?$ have you double checked all of your natural units and scaling? That sort of error magnitude makes it seem like some constants are in the wrong place $\endgroup$
    – whpowell96
    Feb 21 at 20:36

1 Answer 1


Seems that your electron energy of $E=420~\text{MeV}$ messes angle range upon which you calculate cross-section. If you make it smaller by 2 orders of magnitude, like E = 420/10**2, then you will get expected chart (in addition integration will be very fast) :

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.