Stopping Power Calculation

Table of contents

GitHub icon Check the repo.

Introduction

Here we explain a bit of the physics behind this macro.

According to the Particle Data Group, a particle entering a material will have a mean rate of energy loss (also known as stopping power) defined by the (in)famous Bethe-Bloch formula:

-dEdx=Kz2ZA1β212ln2mec2β2γ2TmaxI2-β2 

where dE represents the amount of energy E that is lost by the particle per unit of distance dx traveled inside a material. Since it is a lost, the decrease is represented by a negative sign (-).

What follows after the equal sign is a series of constants and parameters that define the particle and the material. Not all of them, only the ones that have an influence in the loss of energy. In the following table, these variables are introduced (the numbers in parentheses represent the error of the value):

Symbol Definition Units or value
E Incident particle energy MeV[*]
K 4π NA re2 me c2 0.307 075 MeV cm2 / g
z Charge of the incident particle unitless
Z Atomic number of the material unitless
A Atomic mass of the material g / mol
β velocity in terms of the speed of light c unitless
γ Lorentz factor unitless
me c2 Electron mass x c2 0.510 998 918(44) MeV
T Kinetic energy MeV
I Mean excitation energy eV

[*] Electron-volts (eV) are often used as units of energy and mass in physics, using a system of natural units with c and set to 1.

This beautiful formula looks something like this when we put real numbers in it. In this case, the values for muons in copper:

Bethe-Block formula for muons in copper

The code

Translation of the formula

The macros are written in C++. The function stoppingPower() in simulation1.cpp and simulation2.cpp calculates this formula. To translate it to code, first we calculate gamma (the Lorentz factor), which depends on the energy of the particle:

γ = 1+Em

and so we have:

g = 1 + E/m;

We can then play with different particles by changing the value of the mass (m), or launch a particle with different values of energy (E), or combine both to make a study on different particles with different energies.

Next, we calculate beta squared, the speed of the particle:

β2 = 1+1γ2

which translates in the code to:

b2 = 1 - (1/(g*g));

The kinetic energy of the particle is defined by:

T =2mec2β2γ21+2γmem+mem2

since me c2 ≈ 0,511 MeV, and so 2me c2 = 1.022 MeV, we can shorten the expression for T like this:

T = (1.022*b2*g*g) / ( 1 + 2*g*(0.511/m) + pow(0.511/m,2) );

And finally, setting the value of the constant K as double k = 0.307075; (see table), and calling the charge of the particle (z) q, we get to calculate the stopping power:

-dEdx=Kz2ZA1β212ln2mec2β2γ2TmaxI2-β2 

If the mass and energies are given in MeV, the distances in cm, and A in g / mol, the units of the stopping power will be of MeV cm2 / g. We can multiply the stopping power by ρ0, the density of the material in g / cm3, to have it in units of MeV / cm. That way, it will be more straight-forward to calculate the energy loss per unit of travelled distance. Multiplying the stopping power by the travelled distance will automatically give the energy lost by travelling that distance.

The stopping power expression translated into code is:

SP = ((ro*k*q*q*Z)/(A*b2))*(0.5*log( (1.022*b2*g*g*T) / (I*I) ) - b2);

The final subroutine looks like this:

double stoppingPower(double m, double q, double E, double ro, double Z, double A, double I) {

	double k = 0.307075;
	double g, b2, T, SP;

	printf("\n|******** Stopping-Power values for this material ********|\n");

	g = 1 + (E/m);
	printf("\n Gamma  =  %lf", g);

	b2 = 1 - (1/(g*g));
	printf("\n Beta squared  =  %lf", b2);

	T = (1.022*b2*g*g) / (1 + 2*g*(0.511/m) + pow(0.511/m,2) );
	printf("\n Maximum kinetic energy transferred  =  %lf MeV ", T);

	SP = ((ro*k*q*q*Z)/(A*b2))*(0.5*log( (1.022*b2*g*g*T) / (I*I) ) - b2);
	printf("\n Stopping Power  =  %lf MeV/cm \n", SP);

	return SP;
}

Physics and code variables

The simulation is simple: we have a material surrounded by air. At time 0, a particle enters a material of a certain thickness. To track the travelling of the particle in the material, we use steps of a fixed size. Inside a loop, we calculate how much energy the particle loses in each step, and determine if it is stopped inside the material or if it goes out, in which case we repeat the process to calculate the energy loss in air.

Here is a diagram that shows what we have just said:

Diagram of a particle entering a material

The first simulation uses the remaining distance to the exit L:

for(n = 0, L = thickness; L > 0 && E > 0; n++) {
	E = E - SP*step;
	L = thickness - n*step;
}

where n is a counter for the number of steps, and we exit the loop when whether the particle is stopped (E > 0 condition) or it has abandoned the material (L > 0 condition).

The second simulation uses the range:

for(n = 0; range < thickness && E > 0; n++) {
	SP = stoppingPower(m, q, E, ro, Z, A, I, &G, &B, &t);
	E = E - SP*step;
	range += step;
}

Also, the second simulation calculates the stopping power every time the energy changes. The first simulation doesn't because it is intended to be a rough approximation.

Test with real values

Play with the values of energy of the particle and thickness of the material to see what happens. Here is an example of a proton particle entering an aluminium foil. The values can be found in the Wikipedia or in specialized publications, just pay attention to the units and make the necessary conversions:

Parameter Value
Particle values (proton)
Mass (m) 938.272 MeV / c2
Charge (z or q) 1 (in units of e)
Initial kinetic energy (E) 100 MeV[*]
Material values (Al foil)
Density (ρ0) 2.7 g / cm3
Atomic number (Z) 13
Mass number (A) 27
Ionization Potential (I) 0.000006 MeV
Material thickness 1 cm
Step size 0.01 cm[**]
Stopping-Power values for this material
Gamma (γ) 1.106579
Beta squared (β2) 0.183351
Maximum kinetic energy transferred (Tmax) 0.229180 MeV
Stopping Power (- dE / dx) 22.572949 MeV/cm
Result The particle exits the material with an energy of 77.427051 MeV in 100 steps.

[*] You can make a study in energy dependence changing this value systematically.
[**] If more precision is needed, introduce a shorter step size.

If the value of the initial energy is changed, a plot can be made of the dependence of the stopping power with the energy of the incidental particle, or a plot of the energy against the range, like the one below:

Plot of the initial energy against the range

More

More about the code in the README file.

GitHub icon Check the repo.

Twitter icon Share on Twitter or follow me!

Calendar icon Published: 2004-06-01 15:36

Particle Physics Snippets site is coded and maintained by Octopus in Vitro.
You can also find the code on BitBucket and GitLab.