22.1 Reading Simulation Output with Python

Table of Contents

NumPy

In GATE, most numerical data that you want to analyze can be represented as NumPy arrays. You will almost always use NumPy as a first step after reading simulation output, because it gives you fast operations on large collections of numbers.

When you read data, you typically end in one of two situations. Either you directly obtain NumPy arrays, or you use another tool, such as Pandas or uproot, and then convert to NumPy. Once you have arrays, you can compute basic statistics, build histograms, or reshape the data.

Assume you already have three arrays, energy, x, and y, representing for example the energy and position of singles in a detector. You can compute simple quantities in a very compact way:

python

import numpy as np
mean_energy = np.mean(energy)
std_energy = np.std(energy)
max_x = np.max(x)
min_y = np.min(y)

NumPy can also create histograms that you can later plot. For example, to build an energy spectrum:

python

hist, bin_edges = np.histogram(energy, bins=200, range=(0, 1.0))  # MeV

The result hist contains the counts per bin and bin_edges describes the bin boundaries. You can use these arrays directly with Matplotlib to create energy spectra.

For multidimensional data, you often store several columns in a single array of shape (N, M), where N is the number of events and M is the number of variables. Indexing along the second dimension then gives you one quantity:

python

data = np.column_stack([energy, x, y])
energy_column = data[:, 0]
x_column = data[:, 1]

This structure is useful if you read text output, such as CSV files from GATE, with np.loadtxt or np.genfromtxt. For example:

python

data = np.loadtxt("singles.csv", delimiter=",", skiprows=1)
energy = data[:, 0]
time = data[:, 1]

You can then combine masks and conditions to filter the events. For example, to select events in a photopeak window:

python

mask = (energy > 0.45) & (energy < 0.60)  # MeV
photopeak_events = data[mask]

NumPy operations work on entire arrays at once, which is essential when your GATE simulations produce millions of hits or singles. You avoid explicit Python loops, which are slow for large datasets.

You may also want to store intermediate results. NumPy provides simple IO functions for that purpose, for example:

python

np.save("energy.npy", energy)
loaded_energy = np.load("energy.npy")

Or you can store several arrays in one file:

python

np.savez("results.npz", energy=energy, x=x, y=y)
saved = np.load("results.npz")
energy_loaded = saved["energy"]

Pandas

Pandas is useful when you want to treat your GATE output as tabular data, where each column has a name such as energy, time, or crystalID. It adds labeled columns, convenient filtering, and summarizing operations on top of NumPy.

A common workflow is to read CSV or other text-based output from GATE directly into a DataFrame. Assume your simulation wrote a text file with singles, with a header line:

python

import pandas as pd
df = pd.read_csv("singles.csv")

Pandas automatically converts the columns into internal arrays and attaches column names from the header. You can inspect the data:

python

print(df.head())
print(df.columns)
print(df.describe())

This is helpful for quick checks of your simulation output. You can then select subsets of data using conditions. For example, select singles in an energy window and in a specific detector:

python

energy_mask = (df["energy"] > 0.45) & (df["energy"] < 0.60)  # MeV
detector_mask = df["detectorID"] == 5
subset = df[energy_mask & detector_mask]

From this point, you can compute statistics:

python

mean_time = subset["time"].mean()
count_events = subset.shape[0]

If you want to build histograms with Pandas, you can either use its own plotting interface or convert to NumPy. For example:

python

energy_array = df["energy"].to_numpy()

Then you can use NumPy or Matplotlib. You can also group events by detector and compute summary information, which is particularly useful in detector studies:

python

grouped = df.groupby("detectorID")["energy"]
mean_by_detector = grouped.mean()
counts_by_detector = grouped.count()

This gives you detector occupancy maps and energy per detector without writing explicit loops.

Pandas can also write processed data back to disk in CSV or other formats, which is useful if you want to store intermediate analysis stages:

python

subset.to_csv("singles_photopeak_detector5.csv", index=False)

Although Pandas is built on top of NumPy, its tabular interface and named columns often make analysis scripts shorter and easier to read, especially when you have many different variables in your GATE output.

Uproot

Most GATE simulations use ROOT files as the main output format. To analyze these ROOT files with Python, you can use the uproot library. Uproot reads ROOT files in pure Python and can convert ROOT trees into NumPy arrays or Pandas DataFrames without requiring a compiled ROOT installation.

At a basic level, a GATE ROOT file contains one or more trees. Each tree has branches such as EventID, Energy, GlobalTime, or detector identifiers. The main idea is to open the file, select a tree, and then read the branches.

First, install uproot in your Python environment if it is not already installed:

bash

pip install uproot

Then, in a Python script or notebook, you can write:

python

import uproot
file = uproot.open("gate_output.root")
print(file.keys())

This prints the names of the objects in the file. Identify the tree you want, for example "Singles;1" or "Hits;1". Then you access the tree:

python

tree = file["Singles"]
print(tree.keys())

The tree keys are the branch names. To read several branches into NumPy arrays at once, use:

python

arrays = tree.arrays(["energy", "global_time", "crystalID"], library="np")
energy = arrays["energy"]
time = arrays["global_time"]
crystal_id = arrays["crystalID"]

Here energy, time, and crystal_id are standard NumPy arrays, which you can process as described in the NumPy section.

If you prefer Pandas, you can ask uproot to return a DataFrame:

python

import pandas as pd
df = tree.arrays(["energy", "global_time", "crystalID"], library="pd")

You now have a table where columns correspond to ROOT branches. You can filter and group data using the Pandas tools.

GATE ROOT trees may be large. Uproot can read them in chunks to avoid loading the entire dataset into memory. This is especially important for long simulations or detailed phase space output. You can iterate over chunks like this:

python

for arrays in tree.iterate(["energy", "global_time"], step_size=100_000, library="np"):
    energy_chunk = arrays["energy"]
    time_chunk = arrays["global_time"]
    # process the chunk here

You can accumulate histograms or statistics across all chunks without keeping all events in memory at once.

Some branches in ROOT trees can be jagged arrays, that is, they contain a different number of entries per event. Uproot represents these with its own types, but you can convert them to NumPy or manipulate them directly. For most simple GATE outputs such as singles or coincidences you will typically have flat arrays, which work as shown above.

To summarize the typical path for ROOT analysis in Python: open the ROOT file with uproot, inspect the available trees and branches, read the needed branches into NumPy or Pandas, then apply the analysis steps of your choice, for example energy spectra, timing distributions, or detector maps.

Always verify the ROOT branch names and units before analysis. Different GATE examples or versions may use slightly different branch names or units for energy and time. Mismatched assumptions about branch names or units can silently invalidate your analysis.

22.1 Reading Simulation Output with Python

NumPy

Pandas

Uproot

Comments

Where to Move