bethe-bloch.py

"""
K/A = 0.307075 MeV g^-1 cm^2

-dE/dX = K/A z^2 Z / \beta^2 [0.5*ln(2 m_e \beta^2 \gamma^2 T_{max}/I^2) - \beta^2 - \delta(\beta \gama)/2]

z^2 = 1
Z = 1
p = 850 MeV
M = 139.6 MeV
\gamma = sqrt(1+(p/M)^2) = 6.17
\beta = sqrt(1-1/sqr(gamma)) = 0.9868
\beta \gamma = p / M = 6.09
m_e = 0.511 MeV
T_{max} = 2 m_e \beta^2 \gamma^2 / (1 + 2 \gamma m_e / M + (m_e / M)^2)
        \simeq 37.9 MeV
I = 0.00000218 MeV
\rho = 0.0708 g cm^-3
hwp = sqrt(\rho (g cm^-3) Z/A) * 28.816 \simeq  
"""

from math import sqrt, log
sqr = lambda x: x*x

K = 0.307075 # MeV g^-1 cm^2
z = 1 # pion charge
Z_over_A = 0.99212 # hydrogen Z/A
M = 139.6 # MeV is a charged pion mass
m_e = 0.511 # MeV is the electron mass
rho = 0.0708 # g cm^-3 is a hydrogen density
gamma = lambda p: sqrt(1 + sqr(p / M))
beta = lambda p: sqrt(1 - 1 / sqr(gamma(p)))
beta_gamma = lambda p: p / M
T_max = lambda p: 2 * m_e * sqr(beta_gamma(p)) / (1 + 2 * gamma(p) * m_e / M + sqr(m_e / M))
h_omega_p = sqrt(rho * Z_over_A) * 28.816 / 1000 / 1000

dE_over_dX = lambda p: K * sqr(z) * Z_over_A / sqr(beta(p)) * \
    (0.5 * log(2 * m_e * T_max(p) / sqr(h_omega_p)) - sqr(beta(p)) + 1/2)

print beta_gamma(1330), beta(1330)*gamma(1330)
print dE_over_dX(1330)