I have a Geant4 simulation where 662keV photons are going into a detecting volume (Germanium). As expected, when I look at the output of the simulation, these photons undergo Compton scattering and photoelectric absorption.

I have written a python script for analysing this output, where I calculate the scattering angle and plotted them in a histogram:Histogram of Compton scattering angles from a 662keV

This seems to agree very well with plots I've found online. (figure2: DOI:10.1016/j.nima.2019.163228 ) Online plot from a paper.

To round off this piece of work, I wanted to show that my distribution from Geant4 can be inferred from the Klein-Nishina model/the model that the aforementioned paper used.

However, when I plot the differential cross section (which I understood as akin to 'probability'?) vs scattering angle of the KN model, they don’t agree. The distribution shape is totally different.

Klein-Nishina Model

The simulated electrons in the detecting volume aren't free electrons, so I tried adding an Impulse Approximation (Hartree-Fock if that means anything) but it still doesn't line up with the simulation results (although I'm not entirely sure if I did it right – the maths was slightly heavy for me).

How would I go about showing the agreement between theory (Klein-Nishina or other model) and my histogram?

And how did the author of the 'relative probability' graph obtain that data; a model that I don't know of, experimental...?

Any help or push in the right direction is hugely appreciated. :)

EDIT The answer below yields the below python script - although the vertical scale is off.

import numpy as np
import matplotlib.pyplot as plt

c = 299792458.0  # meter / second
E = 662000 *1.602176634 * 10**(-19) #incident energy 
h = 6.62607015 * 10**(-34)  # joule second
hbar = h / (2 * np.pi)
me = 9.1093837015 * 10**(-31)  # kilogram

omega = E / hbar #angular freq
R = hbar * omega / (me * c**2) #variable R in solution

def I(t):#define the function I(theta)
    return -np.cos(t) / R**2 + \
           np.log(1 + R * (1 - np.cos(t))) * (1/R - 2/R**2 - 2/R**3) - \
           1 / (2 * R * (1 + R * (1 - np.cos(t)))**2) + \
           1 / (1 + R * (1 - np.cos(t))) * (-2/R**2 - 1/R**3)

def w_prime(t):#Calc angular_freq for scattered photon
    return omega/(1+(hbar*omega/(me*c**2))*(1-np.cos(t)))

def f(t):#calc f(theta), aka the PDF
    return (1/(I(np.pi)-I(0))) * (w_prime(t)/omega)**2 * (omega/w_prime(t) + w_prime(t)/omega -(np.sin(t)**2))*np.sin(t)

theta_vals = np.linspace(0, np.pi, 1000)

plt.plot(theta_vals, f_vals)
  • $\begingroup$ what formula did you use to go from probability per d-theta to solid angle? $\endgroup$
    – JEB
    Aug 9, 2023 at 23:00
  • $\begingroup$ The bottom plot (with the solid angle) was from a blog on scipy blog, just a few formating changes. scipython.com/blog/the-kleinnishina-formula Its the standard Klein-Nishina model by the looks of it. @JEB $\endgroup$
    – Owain G
    Aug 10, 2023 at 16:32
  • $\begingroup$ Since your plot has a different y axis, why should they look the same? $\endgroup$
    – JEB
    Aug 10, 2023 at 17:01

1 Answer 1


$d\sigma$ needs to be converted to a probability density function $f(\theta)$, then plot $f(\theta)$.

The well-known cross section for Compton scattering is $$ \frac{d\sigma}{d\Omega} =\frac{\alpha^2(\hbar c)^2}{2s} \left(\frac{\omega'}{\omega}\right)^2 \left(\frac{\omega}{\omega'}+\frac{\omega'}{\omega}-\sin^2\theta\right) $$ where $$ \omega'=\frac{\omega}{1+\frac{\hbar\omega}{mc^2}(1-\cos\theta)} $$

We can integrate $d\sigma$ to obtain a cumulative distribution function. Let $I(\theta)$ be the following integral of $d\sigma$. (The $\sin\theta$ is due to $d\Omega=\sin\theta\,d\theta\,d\phi$) \begin{equation*} I(\theta)= \int \left(\frac{\omega'}{\omega}\right)^2 \left(\frac{\omega}{\omega'}+\frac{\omega'}{\omega}-\sin^2\theta\right) \sin\theta\,d\theta \end{equation*}

The solution is \begin{multline*} I(\theta)=-\frac{\cos\theta}{R^2} +\log\big(1+R(1-\cos\theta)\big)\left(\frac{1}{R}-\frac{2}{R^2}-\frac{2}{R^3}\right) \\ {}-\frac{1}{2R\big(1+R(1-\cos\theta)\big)^2} +\frac{1}{1+R(1-\cos\theta)}\left(-\frac{2}{R^2}-\frac{1}{R^3}\right) \end{multline*} where \begin{equation*} R=\frac{\hbar\omega}{mc^2} \end{equation*}

The cumulative distribution function is \begin{equation*} F(\theta)=\frac{I(\theta)-I(0)}{I(\pi)-I(0)}, \quad 0\le\theta\le\pi \end{equation*}

The probability of observing scattered photons in the interval $\theta_1$ to $\theta_2$ is \begin{equation*} P(\theta_1\le\theta\le\theta_2)=F(\theta_2)-F(\theta_1) \end{equation*}

Differentiate $F(\theta)$ to obtain normalized probability density function $f(\theta)$. $$ f(\theta)=\frac{dF(\theta)}{d\theta} =\frac{1}{I(\pi)-I(0)} \left(\frac{\omega'}{\omega}\right)^2 \left(\frac{\omega}{\omega'}+\frac{\omega'}{\omega}-\sin^2\theta\right) \sin\theta $$

Note that if we had carried through the $\alpha^2(\hbar c)^2/2s$ in $I(\theta)$, it would have canceled out in $F(\theta)$.

  • $\begingroup$ I don't really understand how this answer can be correct. I(theta) is everywhere negative, but what is being integrated is positive definite $\endgroup$
    – Knox Long
    Sep 7, 2023 at 14:49

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.