26.9. Plotting the Depth-Dose Curve

Table of Contents

Understanding the Depth‑Dose Curve

In a proton beam in water simulation you usually end with a set of energy depositions in thin slices of the water phantom. The depth‑dose curve is a plot of dose (or deposited energy) as a function of depth along the beam axis. It shows the entrance region, the build‑up, the Bragg peak and the rapid fall‑off beyond the peak.

Your Geant4 application already records energy deposition per slice and you have written this information to a file, for example as an ntuple using G4AnalysisManager. Plotting the depth‑dose curve is mostly an analysis task that happens after the simulation has finished.

From Energy Deposition to Dose

The simulation typically stores, for each event and for each slice, the total energy deposited. To convert this to dose you need the mass of the slice. Consider a slice that is a box of area $A$ and thickness $\Delta z$ with material density $\rho$. The mass is
$$
m = \rho \, A \, \Delta z .
$$

If $E_{\text{dep}}$ is the energy deposited in that slice, the dose $D$ in gray is
$$
D = \frac{E_{\text{dep}}}{m} .
$$

In Geant4 units energy is in joules if you use SI units, but in practice you usually store MeV. You must convert to joules before dividing by the mass in kilograms.

Dose calculation rule
If $E_{\text{dep}}$ is stored in MeV, the slice volume is $V$ in $\text{cm}^3$, and the density is $\rho$ in $\text{g}/\text{cm}^3$, then the dose in gray is
$$
D \,[\text{Gy}] = \frac{E_{\text{dep}} \,[\text{MeV}] \times 1.602\times10^{-13} \,[\text{J/MeV}]}{\rho \,[\text{g/cm}^3] \times V \,[\text{cm}^3] \times 10^{-3} \,[\text{kg/g}] } .
$$

Often you are interested in the shape of the curve, not its absolute normalization. In that case it is common to normalize the dose values so that the maximum is 1. The normalized depth‑dose curve then becomes

$$
D_{\text{norm}}(z) = \frac{D(z)}{D_{\max}} .
$$

This is also convenient when comparing different beam energies or physics settings.

Binning the Data in Depth

To plot dose as a function of depth you need a depth coordinate for each slice. If you created equal‑thickness slices along the $z$ axis with thickness $\Delta z$, you can assign a depth to each slice by its center position.

If the $i$‑th slice covers $z \in [z_i, z_i + \Delta z]$ and the beam enters from the negative $z$ side, a simple choice is the center

$$
z_{\text{center},i} = z_i + \frac{\Delta z}{2} .
$$

For a phantom that starts at $z = 0$ and extends to $z = L$, with $N$ slices of equal thickness, you have

$$
\Delta z = \frac{L}{N}, \quad
z_{\text{center},i} = \left(i + \frac{1}{2}\right) \Delta z, \quad i = 0,1,\dots,N-1 .
$$

In most analyses you sum the energy deposition over all events for each slice, then divide by the number of primary protons to obtain the mean deposited energy per primary per slice. If you want dose per primary, apply the mass division and any unit conversions after that averaging.

A typical data table for plotting might look like:

Slice index $i$	Center depth $z_{\text{center},i}$ (mm)	Mean $E_{\text{dep}}$ per primary (MeV)	Dose (Gy / primary)
0	1.0	0.0012	$D_0$
1	3.0	0.0015	$D_1$
…	…	…	…
$N-1$	$L - \Delta z/2$	0.0000	$D_{N-1}$

You will plot either the dose or the normalized dose against the center depth.

Plotting with External Tools

Although Geant4 can create histograms directly, depth‑dose curves are usually inspected with external analysis tools such as ROOT or a plotting library in Python. The simulation writes one row per slice per event or one row per slice after accumulation, and the analysis script does the rest.

If you wrote a ROOT file with an ntuple that contains per‑slice quantities, for example sliceID, z, and Edep, you can create the depth‑dose array by looping over entries, summing Edep per sliceID, then converting to dose. Finally you create a graph of D(z) or D_norm(z). The important point is to keep all relevant geometry parameters like slice thickness and density consistent between the simulation and the analysis so that the dose calculation is correct.

If you instead wrote a simple text file or CSV, you can read it with Python or another environment, compute dose using the formulas above, and create the curve with any standard plotting library.

Reading and Interpreting the Curve

Once the depth‑dose curve is plotted, you can visually identify how the proton beam behaves in water. Near the entrance, the dose is relatively low and increases gradually as the protons slow down. Near the end of their range you see the Bragg peak, a sharp maximum where most of the energy is deposited. Beyond that peak the dose drops quickly to almost zero, showing one of the key advantages of proton therapy.

For quantitative work you may want to extract specific points from the curve, such as:

The depth at which the dose reaches its maximum, the Bragg peak position.
The distal fall‑off, for example the depths where the dose drops to 80 percent and 20 percent of the maximum value.
The entrance dose compared to the peak dose.

These quantities are computed by scanning the discrete depth‑dose points from the plot and applying criteria such as a percentage of $D_{\max}$.

If you run simulations with different beam energies, you will obtain different curves. Higher energy protons have a deeper Bragg peak. By plotting several curves on the same graph with clear labels you can compare how the range and peak shape change with energy, physics list choice, or phantom material.