KAHIBARO
22.3. Position Distributions
Table of Contents
Interaction coordinates
In GATE, every interaction of a particle with matter can be recorded along with its spatial coordinates. For position distributions, you typically use outputs that contain per-interaction or per-event positions, such as hits, singles, coincidences, or phase-space data. In this chapter we focus on how to handle these positions in Python so you can understand where interactions occur inside your simulated system.
When you read GATE output with Python, for example using uproot and NumPy or Pandas, interaction coordinates are usually stored as three separate columns or branches, such as x, y, and z, or sometimes as globalPosX, globalPosY, globalPosZ. The units follow the GATE convention, most often millimeters. It is essential to verify the unit of the stored positions in the documentation or by checking your simulation script. A common pattern is that you wrote something like my_actor.spacing = [1 mm, 1 mm, 1 * mm] for images or that geometry dimensions were given in millimeters. In analysis you must treat the numeric values accordingly, for example converting from millimeters to centimeters when plotting dose as a function of depth.
You will typically load the ROOT file with uproot, extract the position arrays, and convert them into NumPy arrays. From there, you can compute distributions of positions. One frequent analysis is the 1D distribution of interaction depth along a chosen axis. For instance, in a water phantom illuminated by a proton beam, you might look at the histogram of z positions where energy is deposited. With Python this becomes a simple call to numpy.histogram on your z array, using an appropriate binning. For PET or SPECT, you may do something similar with the radial coordinate, for example $r = \sqrt{x^2 + y^2}$, to see how often photons interact at different distances from the center.
To analyze interaction coordinates effectively, you should think in terms of the geometry coordinate system of your simulation. The world coordinate system is usually centered at the origin, with the $z$ axis along the scanner axis or beam axis, but this depends on how you created the geometry. When you create position distributions, you should relate each axis to a meaningful direction, such as radial, axial, or depth.
You can also transform coordinates to more convenient systems for analysis. For example, for cylindrical PET scanners, it is often useful to convert from $(x, y)$ to polar coordinates $(r, \phi)$, where
$$r = \sqrt{x^2 + y^2},$$
and
$$\phi = \arctan2(y, x).$$
With NumPy, you can compute these for all events in one operation. This lets you study how interactions are distributed around the ring or how uniform the system is in angle.
In many analyses, you will correlate position with other quantities. One common example is a 2D histogram of position versus energy. For example, for a scintillator block, you might plot a 2D density of $(x, y)$ locations colored by the mean deposited energy to see if the energy response is uniform over the surface of the detector. Another is time versus position, for example the interaction time of photons as a function of depth to investigate time-of-flight effects.
For voxelized dose or energy deposition, you often do not work with lists of individual interaction coordinates, but with pre-binned 3D images produced by actors that already summarized energy or dose per voxel. Even in that case, understanding interaction coordinates helps you interpret where each voxel center is located in space. The index $(i, j, k)$ in an array corresponds to a physical coordinate
$$x = x_0 + (i + 0.5) \cdot \Delta x,$$
$$y = y_0 + (j + 0.5) \cdot \Delta y,$$
$$z = z_0 + (k + 0.5) \cdot \Delta z,$$
where $(x_0, y_0, z_0)$ is the origin and $(\Delta x, \Delta y, \Delta z)$ are the voxel sizes. These relationships are often provided by the metadata that accompanies MHD or NIfTI images.
Always ensure that:
- You know the coordinate system (origin and axis directions) used in your simulation.
- You confirm the units of position values before plotting or computing distances.
- You consistently convert between index space $(i, j, k)$ and physical coordinates $(x, y, z)$ when working with voxel images.
When making 1D, 2D, or 3D histograms of interaction coordinates in Python, take care in the choice of bin size and range. Too few bins will hide details, while too many bins will create noisy distributions, especially for simulations with a limited number of events. You can estimate the bin size from typical physical dimensions, for example a bin width similar to crystal size for detector maps or to voxel spacing for dose distributions.
Finally, position distributions can also be used for quality checks of your simulation. If you know that all interactions should happen in a particular volume, such as a phantom or detector, and you observe a significant number of events at positions outside that region, this indicates a problem in geometry or actor configuration that you should correct before proceeding.
Detector maps
Detector maps are spatial summaries of interaction positions, usually projected into a 2D plane associated with a detector surface or volume. Their purpose is to show how often and where inside a detector interactions or detected events occur. Detector maps are especially important for imaging detectors, such as PET crystals, gamma camera crystals, or CT detector arrays.
In Python, you typically start from a list of events with associated detector identifiers and possibly positions. Depending on how your simulation is configured, you might have one or both of these: direct $(x, y, z)$ positions of hits or singles, and integer IDs such as crystalID, moduleID, or pixelID that encode which detector element was involved. You can build detector maps directly from coordinates using 2D histograms, or from IDs using integer-based aggregation.
A simple detector map based on positions uses a 2D histogram of $x$ versus $y$ at a fixed $z$ plane, corresponding to the detector surface. With NumPy and Matplotlib, you can compute a 2D histogram of arrays x and y, then display it as an image. The color of each pixel in the image represents the number of interactions in that spatial bin. By choosing the bin size to match the physical segmentation of your detector, such as the crystal width or pixel pitch, you can see uniformity, dead regions, and edge effects.
When your data are better represented with discrete detector IDs, a detector map becomes a 2D array indexed by these IDs. For example, in a gamma camera, you may have a rectangular array of pixels with indices $(i, j)$. You can create an array with shape (n_rows, n_cols) and increment its entries based on the pixel indices from your data. Similarly, for PET systems, you may group events by crystal index in one ring module and create a 2D map of radial versus axial crystal index.
Detector maps can include not only counts but also derived quantities. For each detector element or spatial bin, you might record the mean energy, standard deviation of energy, or other metrics. Python tools like Pandas make this straightforward. For example, you can group your table of events by crystalID and compute count, mean_energy, and std_energy for each group, then reshape these values into 2D arrays that match the physical layout of the detector.
A useful pattern is to store the detector layout in a small table that maps from detector IDs to geometric coordinates. For example, you might maintain a table with columns crystalID, x_center, and y_center. Once loaded in Python, you can merge this information with your event data, and then create scatter plots or interpolated maps where each crystal center is plotted at its physical location with color representing counts or mean energy. This is valuable when the detector layout is not a simple regular grid, such as ring segments in PET or complicated block detectors.
The table below summarizes common ways to build detector maps from GATE output using Python.
| Data available | Typical Python approach | Detector map type |
|---|---|---|
| $(x, y)$ positions | numpy.histogram2d(x, y, bins=...) | 2D spatial histogram in mm |
| Pixel indices $(i, j)$ | 2D NumPy array counts[i, j] += 1 | Pixel occupancy map |
| Crystal IDs | Groupby with Pandas, then reshape | Crystal occupancy or energy maps |
| Crystal ID + coordinates | Merge on ID, then scatter with Matplotlib | Non regular detector geometry map |
Detector maps are especially powerful for detector validation. For a uniform irradiation, you expect a uniform detector map, possibly with edge effects. If your map shows unexpected hot or cold spots, repeating patterns, or missing areas, this may indicate geometry problems, misconfigured sensitive volumes, or issues in the digitizer chain that discards some events.
For PET and SPECT, you can also extend the concept of detector maps to show coincidences or projection counts. For example, in PET you might build a 2D histogram of detected coincidences as a function of two crystal indices, which resembles a sinogram. For SPECT, you can sum detected counts over energies and time to create a 2D projection image, which is itself a detector map in the detector coordinate system.
In all cases, when you produce detector maps in Python you must ensure that the mapping from IDs or indices to physical positions is correct and consistent with the geometry defined in your GATE simulation. Store this mapping explicitly in your project, for example in a configuration file or a small helper module, so that future analyses and collaborators can reproduce your detector maps without confusion.
Detector maps also form an intermediate step to more advanced performance analysis. By combining maps for different conditions, such as different energy windows or time intervals, you can study position dependent detector efficiency, spatial variations in energy resolution, or timing performance across modules. These analyses are central to evaluating and improving the design of simulated detectors and to validating the realism of your GATE models against experimental data.
Views: 12