26.7. Calculating Depth Dose

Table of Contents

Understanding Depth Dose

Depth dose describes how much energy a beam deposits as a function of depth inside a material. For a proton beam in water, this curve shows how the dose rises, peaks at the Bragg peak, and then falls off sharply. In the simulation, you already divided the water phantom into slices and recorded the energy deposited in each slice. Calculating depth dose means turning those energy depositions into a quantity that can be interpreted as dose per slice and plotting it versus depth.

In Geant4, this step is mainly an analysis task. The geometry, tracking, and scoring of energy deposition are already in place, and you now combine those results into a depth dependent observable.

From Energy Deposition to Dose

Dose $D$ is defined as the absorbed energy per unit mass. In SI units it is measured in gray, where $1\ \text{Gy} = 1\ \text{J/kg}$.

For each phantom slice, you have the total energy deposition from all steps in that slice, usually in units like MeV. You can convert this to dose by dividing by the mass of the slice and applying the unit conversion between MeV and joules.

A typical chain is:

For each slice, accumulate the total deposited energy $E_{\text{dep}}$ during an event (or during the whole run).
Compute the slice mass $m$ from its volume and density.
Convert $E_{\text{dep}}$ to joules.
Compute dose $D = E_{\text{dep}} / m$.

Key relation for depth dose
The absorbed dose in a slice is
$$
D = \frac{E_{\text{dep}}}{m}
$$
where:

$D$ is dose in gray (Gy),
$E_{\text{dep}}$ is the deposited energy in joules (J),
$m$ is the mass of the slice in kilograms (kg).
To convert energy from MeV to joules:
$$
1\ \text{MeV} = 1.602 \times 10^{-13}\ \text{J}.
$$

In many practical proton depth dose simulations, you are interested in the shape of the curve rather than the absolute dose. In that case, you can simply normalize the curve to its maximum or to the entrance value and treat the vertical axis as relative dose. The calculation process in Geant4 is the same, but you can skip the physical unit conversion step if you only want relative values.

Using Slice Geometry to Compute Mass

To compute dose, you need the mass of each slice of the water phantom. If every slice has the same thickness and cross section, they all have the same volume and therefore the same mass.

Let the slice dimensions be length $L_x$, width $L_y$, and thickness $L_z$:

The slice volume is $V = L_x L_y L_z$.
For water, the density is approximately $1.0\ \text{g/cm}^3$ in human terms. In Geant4 units, density is typically specified in $\text{g}/\text{cm}^3$ as well, but the underlying system is based on $\text{kg}/\text{m}^3$.

You can obtain the slice volume directly from your geometry. For a simple box slice:

If you explicitly created each slice as a separate solid, its dimensions are known when you build the geometry.
If you used replica or parameterized volumes, the slice thickness is specified in the replication or parameterization, and the common cross section is the phantom cross section.

In code, you can compute the mass value once and reuse it for all slices. For water:

Density: $\rho \approx 1.0\ \text{g/cm}^3 = 1000\ \text{kg/m}^3$.
Volume: for example, $L_x = 10\ \text{cm}$, $L_y = 10\ \text{cm}$, $L_z = 1\ \text{mm}$, converted to consistent units.

You can either perform the calculation numerically outside Geant4, based on known dimensions, or, more flexibly, retrieve the logical volume of one slice and query its solid for the volume and its material for the density. This latter method keeps your analysis independent of hard coded geometry values.

Once you have the volume $V$ (in $\text{mm}^3$ or $\text{cm}^3$) and density $\rho$, you compute the mass as $m = \rho V$ and convert to kilograms if needed for SI dose units.

Event by Event vs Run Totals

There are two common ways to compute depth dose in a simulation, and both are useful for different purposes.

In one approach, you compute the dose per slice per event. For each event, you accumulate the total deposited energy in each slice during that event, and at the end of the event you convert to dose and fill a two dimensional data structure such as an ntuple with columns for depth and dose. This allows you to later analyze event by event fluctuations and statistical properties of the dose at each depth.

In another approach, you accumulate energy deposition over the entire run before converting to dose. In this method, you sum the deposited energy in each slice over all events in a single accumulator per slice. At the end of the run, you divide each accumulator by the number of primary particles (or the total number of events) to obtain the average energy deposited per primary, convert to dose per primary, and then if desired apply any normalization. This method gives you directly the mean depth dose curve over the whole run.

In both cases, you need to ensure that the per slice accumulators are properly initialized at the start of each run or each event, and that they are used consistently. Within Geant4, RunAction is a natural place to initialize and finalize run level accumulators, while EventAction is the place where event level sums are reset and written out at the end of each event.

Depth Binning and Slice Index

To build a depth dose curve, you need to associate each energy deposition with a depth. You already have a phantom divided into slices, so the slice index directly corresponds to a depth bin.

In practice, each time a step occurs, you identify the slice in which the step happened. This is typically done by:

Assigning unique copy numbers to each slice when you place or replicate them, and
Fetching the copy number from the touchable handle in SteppingAction or SensitiveDetector.

The copy number then acts as a depth bin index. You map copy number $i$ to a depth position, for example the center of the $i$th slice. If each slice has a thickness $\Delta z$, the depth coordinate for slice $i$ can be defined as:

$$
z_i = z_{\text{start}} + \left(i + \frac{1}{2}\right) \Delta z
$$

where $z_{\text{start}}$ is the position of the front face of the phantom along the beam direction.

You may store both the bin index and its corresponding depth value in your analysis output, but for plotting the depth dose curve later it is usually enough to know the bin order and thickness, and reconstruct the exact depth axis in your analysis tool.

Normalizing the Depth Dose Curve

Once you have dose or average deposited energy per slice, you often want to normalize the curve. Normalization makes it easier to compare different simulations or to compare with experimental or reference data.

Common normalization choices include:

Normalizing to unity at the maximum: divide all dose values $D_i$ by the maximum dose $D_{\text{max}}$ across all slices, so that the Bragg peak has value 1.
Normalizing to the entrance dose: divide all $D_i$ by the dose in the first slice, to show the relative build up and the peak compared to the entrance.
Normalizing by the number of primary particles: useful if you record absolute energy deposition without already dividing by the number of events.

In code, normalization is typically done at the end of the run, when all per slice accumulators are finalized. You scan for the maximum, or select the entrance slice value, and then rescale all per slice values. You then store or output the normalized values as relative dose versus depth.

Important depth dose normalization practices

Always note which normalization you use, entrance dose or peak dose.
When comparing with other simulations or measurements, use the same normalization convention.
If you use relative dose, state that absolute dose units are not preserved.

Handling Statistical Uncertainties

Since a Geant4 simulation is a Monte Carlo calculation, the dose in each slice is subject to statistical fluctuations. The uncertainty decreases as you increase the number of primary particles.

To quantify the uncertainty, you can compute the mean and variance of the deposited energy or dose in each slice over events. One simple method is to:

For each slice, keep two accumulators: the sum of deposited energy per event, and the sum of the squares of the deposited energy per event.
At the end of the run, compute the mean $\bar{E}$ and variance $\sigma^2$ using:

$$
\bar{E} = \frac{1}{N} \sum_{k=1}^{N} E_k
$$

$$
\sigma^2 = \frac{1}{N} \sum_{k=1}^{N} E_k^2 - \bar{E}^2
$$

where $N$ is the number of events and $E_k$ is the energy deposition in that slice in event $k$.

You can propagate the uncertainty from deposited energy to dose by dividing both the mean and the standard deviation by the same mass factor, since dose is proportional to energy deposition for a fixed mass. This provides an error bar for each point in your depth dose curve, which is important when validating the simulation or comparing different configurations.

Exporting Depth Dose for Plotting

The final step in calculating depth dose for your proton beam in water is exporting the depth dependent values to an output format suitable for plotting and comparison.

If you use the Geant4 analysis system, you typically create an ntuple or histogram with one axis for depth and one for dose or energy deposition. For depth dose curves, an ntuple with columns like "depth", "dose" or "relativeDose", and optionally "error" is often convenient. Each slice contributes one row, written once per run during RunAction at the end of the run.

Outside Geant4, you can read this data with ROOT or another plotting tool and construct a graph of dose vs depth. If you used slice indices only, you reconstruct the depth axis based on the slice thickness and starting position. If you stored explicit depth values, you can plot directly.

By structuring the calculation in this way collecting deposited energy per slice, converting to dose, normalizing, handling uncertainties, and exporting per slice values you obtain a clear and reproducible depth dose curve that reflects the underlying physics of proton interactions in water and that serves as the basis for observing the Bragg peak in the next step.