A hands-on intro to Monte Carlo radiation simulation and analysis
+
+
+
+
+
+
+
+
Author
+
+
John Murphy and Steven Carlson
+
+
+
+
+
+
+
+
+
+
+
+
+
+
What this is
+
A four-lab series that takes you from a fresh Linux install to a working PET-scanner simulation, using Geant4 for the physics and Python (pandas, matplotlib) for analysis.
+
The labs are designed to be done in order. Each one builds on the environment, code, and analysis idioms of the previous one.
+
+
+
Audience
+
Upper-level undergraduates or early graduate students with:
+
+
Comfort with the Linux command line
+
Some C++ literacy (you don’t need to write much; you’ll mostly read and patch)
+
Some Python literacy (you’ll write short analysis scripts)
+
+
No prior Geant4 or medical-physics knowledge is assumed.
Lab 3 — Acquiring and Analysing Simulation Data Patch B1’s C++ to write per-event energy data to CSV via G4AnalysisManager, then load it in Python with pandas and produce histograms, statistics, and per-particle comparisons.
+
Lab 4 — A PET Scanner Apply the same workflow to Geant4’s B3a (Positron Emission Tomography) example. Identify the 511 keV photopeak in the energy spectrum, find true coincidences, and reconstruct lines of response.
+
+
+
+
What you’ll build
+
By the end of the series you will have:
+
+
A reproducible Geant4 + Python development environment
+
Modified Geant4 applications that emit clean CSV data
+
Python scripts that compute statistics and produce publication-style plots
+
A working understanding of the simulate → save → analyse → visualise loop that underlies most computational physics
In this lab, you will set up the development environment that the rest of the course depends on. By the end of this session, you will be able to:
+
+
Install Ubuntu Linux in a virtual machine and configure it for productive use
+
Install the system libraries required to build Geant4 with Qt visualization support
+
Compile Geant4 from source using CMake
+
Configure shell environment variables so applications can locate the Geant4 installation
+
Verify the installation by inspecting the directory layout and environment
+
+
+
+
Introduction
+
In this lab, you will install a Linux operating system in a virtual machine and build Geant4 from source with Qt (visualization) support. Geant4 is a toolkit for simulating ionizing particle interactions with matter — the foundation we will use throughout the rest of the lab series to model detectors, beams, and medical imaging systems.
+
+
+
Prerequisites
+
+
Computer with virtualization support enabled in BIOS/UEFI
At least 50 GB of free disk space and 8 GB of RAM available to the VM
+
+
+
+
Pre-Lab Questions
+
+
What is a virtual machine, and why are we using one for this course rather than installing Linux on the host machine directly?
+
What is the difference between a source distribution of a software package and a binary (pre-compiled) distribution?
+
What does CMake do? Why do we configure a build with CMake before running make?
+
What is an environment variable, and why does Geant4 ship a script (geant4.sh) for setting them?
+
+
+
+
Step 1: Setting Up the Virtual Machine
+
+
+
+VirtualBox New VM
+
+
+
+
Open VirtualBox or VMware.
+
Click “New” to create a virtual machine.
+
Name your virtual machine and select Linux type/version (Ubuntu).
+
+
Do not use the default VM location (create a new accessible folder instead).
+
+
Copy the ISO file into your chosen accessible folder (e.g., Downloads).
+
Select the ISO image file when prompted by VirtualBox.
+
Allocate memory (RAM): at least 8 GB.
+
Allocate CPU cores (maximum 60% of available cores). Check core count in Windows via Settings > About Your PC.
+
Create a virtual hard disk (at least 50 GB recommended).
+
Start the virtual machine and boot from the ISO.
+
Select “Try or Install Ubuntu” to launch the live session. This boots Ubuntu into a live session, allowing you to test it without affecting the virtual disk.
+
From the Linux desktop, click “Install Ubuntu” and follow installation prompts to create a permanent installation.
+
After installation, remove the ISO:
+
+
+
Go to the top menu bar and select Devices > Optical Drives > Remove disk. (Note: if its greyed out, then it has already been removed automatically).
+
Restart the virtual machine.
+
+
+
Log in to your new Ubuntu installation and switch to X11 session at the login screen (click the gear icon).
+
+
+
+
Step 2: Enable Shared Clipboard and Fullscreen (Optional)
+
+
Open Terminal in Linux (Ctrl+Alt+T). Update and install essential packages:
echo"/usr/local/lib"|sudo tee \
+ /etc/ld.so.conf.d/geant4.conf
+sudo ldconfig
+
Set environment variables for Geant4. The Geant4 installer ships a shell script that exports G4INSTALL, G4LEDATA, and the other variables Geant4 applications need at runtime:
+
source /usr/local/bin/geant4.sh
+
You will need to source this script in every new terminal where you build or run Geant4 code. (You can add source /usr/local/bin/geant4.sh to your ~/.bashrc to make it automatic.)
+
+
+
Verification and Analysis
+
+
Verify that the environment is set correctly. In the same terminal where you sourced geant4.sh, run:
+
env|grep G4
+
You should see several lines of output that start with G4... — for example G4INSTALL=/usr/local, G4LEDATAPATH=..., and so on. This confirms the environment variables were exported.
+
If G4INSTALL is unset, the source step did not run successfully; re-run it before continuing.
+
Navigate to your installation directory and explore the structure:
+
cd /usr/local
+ls-R|less
+
Can you find:
+
+
The lib directory containing the libraries?
+
The include directory with header files?
+
The share directory with examples and data files?
+
+
+
+
+
Troubleshooting
+
+
If CMake fails, check that all required packages are installed
+
If compilation fails with memory errors, reduce the number of parallel jobs (e.g., make -j2 instead of make -j$(nproc))
+
If environment variables are not set, ensure you have sourced geant4.sh in your current terminal
+
+
+
+
Verification Checklist
+
Before completing this lab, verify you can:
+
+
+
+
+
+
+
+
+
+
+
Conclusion
+
You have successfully set up a Linux virtual machine and compiled Geant4 from source with Qt visualization support. In Lab 2 you will use this environment to build and run your first particle-physics simulation.
Lab 2: First Simulation - Building and Running Geant4 Example B1
+
+
+
+
+
+
+
+
Author
+
+
John Murphy and Steven Carlson
+
+
+
+
+
Published
+
+
November 11, 2025
+
+
+
+
+
+
+
+
+
+
+
+
+
Learning Objectives
+
In this lab, you will compile and run your first Geant4 simulation. You will learn the difference between interactive and batch mode and how to control the simulation parameters. By the end of this session, you will be able to:
+
+
Locate and copy Geant4 example applications from the source distribution
+
Compile a Geant4 example application using CMake
+
Run a simulation and visualize the results in a graphical user interface (GUI)
+
Use the Geant4 command-line interface to control the simulation
+
Create and use macro files to run simulations in batch mode
+
Modify simulation parameters, such as particle type, energy, and the number of events
+
+
+
+
Introduction
+
Now that you have a working Geant4 installation, it’s time to use it. Geant4 is a toolkit, which means we use it to build our own specific applications. The Geant4 distribution includes many examples that demonstrate how to do this. We will start with Example B1, a simple application that simulates a particle beam hitting a target made of different materials and measures the energy deposited.
+
+
+
+Geant4 Example B1 Geometry
+
+
+
You will learn to run this application in two ways:
+
+
Interactive Mode: This mode opens a GUI, allowing you to see the detector geometry and particle tracks. It’s excellent for debugging and understanding the simulation visually.
+
Batch Mode: This mode runs without a GUI, executing a series of commands from a text file called a “macro”. This is how you run simulations to collect large amounts of data for analysis.
+
+
Mastering both modes is key to the scientific computing workflow: explore and debug visually, then collect data efficiently.
+
+
+
Prerequisites
+
Before you begin, ensure you have the following:
+
+
Completed Lab 1 (Geant4 installed and environment configured)
+
Geant4 source code available on your system
+
CMake and C++ compiler (already installed in Lab 1)
+
Qt5 support enabled (for GUI visualization)
+
+
+
+
Pre-Lab Questions
+
+
In the context of a simulation, what is an “event”? What is a “run”?
+
What is the purpose of a macro file in Geant4?
+
Based on the Example B1 README, what is the default particle and energy being simulated?
+
What is the difference between compiling Geant4 itself and compiling a Geant4 application?
+
+
+
+
Step 1: Locate and Copy Example B1
+
The Example B1 source code (and all other examples) are included in the Geant4 source code that you downloaded in Lab 1.
+
+
First, verify where your Geant4 source code is located. If you followed Lab 1, it should be in your home directory:
+
ls ~/geant4-v11.3.2/examples/basic/B1
+
Create a dedicated directory for lab examples and navigate into it:
+
mkdir-p ~/g4-labs/lab2
+cd ~/g4-labs/lab2
+
Copy the source code for Example B1 into your lab directory:
+
cp-r ~/geant4-v11.3.2/examples/basic/B1 .
+
Navigate into the B1 example directory:
+
cd B1
+
Explore the contents of the directory:
+
ls-l
+
You should see several directories and files:
+
+
src/ - Source code files (.cc)
+
include/ - Header files (.hh)
+
exampleB1.cc - Main program file
+
CMakeLists.txt - Build configuration
+
*.mac - Macro files for different run configurations
+
README - Documentation for this example
+
+
+
+
+
+B1 Directory Structure
+
+
+
+
+
Step 2: Compile the B1 Application
+
Just like we compiled Geant4 itself in Lab 1, we need to compile the B1 application. We will use the same out-of-source build technique with CMake.
+
+
First, ensure your Geant4 environment is set up. Open a new terminal and run:
+
source /usr/local/bin/geant4.sh
+
Verify the environment variables are set:
+
env|grep G4
+
Create a build directory for the B1 application inside the B1 folder:
+
mkdir build
+cd build
+
Run cmake to configure the build. The .. tells cmake to look for the source code in the parent directory:
+
cmake ..
+
Because you sourced the geant4.sh script, cmake automatically knows where to find your Geant4 installation.
+
+
+
+
+CMake Configuration Output
+
+
+
+
Compile the application using make:
+
make-j$(nproc)
+
This uses all available CPU cores to speed up compilation.
+
After compilation finishes, verify the executable was created:
+
ls-lh exampleB1
+
You should see the exampleB1 executable file.
+
+
+
+
+Successful Build Output
+
+
+
+
+
Step 3: Run in Interactive Mode (GUI)
+
Now let’s run the simulation with the graphical user interface to visualize what’s happening.
+
+
From the B1/build directory, run the executable with no arguments:
+
./exampleB1
+
A Qt GUI window should appear showing:
+
+
Left pane: 3D viewer with the detector geometry
+
Right pane: Control panel with viewing options
+
Bottom pane: Command-line interface with Idle> prompt
+
+
+
+
+
+B1 Interactive GUI Window
+
+
+
+
The simulation is now waiting for your commands. Let’s fire 10 particles. At the Idle> prompt, type:
+
/run/beamOn 10
+
You will see green lines appear in the detector. These are the tracks of the particles (6 MeV gamma rays by default) and their secondary particles as they travel through the water, tissue, and bone materials.
+
+
+
+
+Particle Tracks Visualization
+
+
+
+
Experiment with the visualization:
+
+
Use your mouse to rotate the view (left-click and drag)
+
Zoom in/out (scroll wheel)
+
Pan the view (middle-click and drag, or shift + left-click)
+
Try different viewing angles to understand the geometry
+
+
Run more events to see different particle behaviors:
+
/run/beamOn 50
+
+
+
+
+Multiple Events Display
+
+
+
+
To exit the application, type exit at the prompt:
+
exit
+
+
+
+
Step 4: Understanding the Simulation Output
+
+
Look at the terminal output from your interactive run. You should see information including:
+
+
Run summary
+
Mean energy deposited in different volumes
+
Dose calculations
+
+
+
+
+
+Terminal Output Example
+
+
+
+
The output shows:
+
+
Absorber: Energy deposited in the main target material
+
Gap: Energy deposited in the gap regions
+
Dose: Radiation dose in the sensitive detector
+
+
+
+
+
Step 5: Run in Batch Mode (Macros)
+
Running simulations by typing commands one by one is tedious. For longer runs and reproducible simulations, we use macro files.
+
+
Go back to the main B1 directory (one level up from build):
+
cd ..
+
View the contents of the provided macro file run1.mac:
+
cat run1.mac
+
You will see a list of Geant4 commands, including /run/beamOn 10.
+
+
+
+
+Macro File Contents
+
+
+
+
Now, run the simulation using the macro file:
+
cd build
+./exampleB1 ../run1.mac
+
The simulation will run immediately, printing information to the terminal, and then exit. You will see output summarizing the run, including the mean energy and dose deposited in the sensitive volumes. No GUI will appear—this is batch mode.
+
Let’s examine another macro file that shows more options:
+
cat ../exampleB1.in
+
Notice the lines controlling the particle gun:
+
+
/gun/particle gamma - Sets particle type
+
/gun/energy 6 MeV - Sets particle energy
+
+
+
+
+
Step 6: Modifying the Simulation
+
The real power of Geant4 is the ability to easily change simulation parameters using UI commands in macro files.
+
+
Create a new macro file for your custom simulation. In the B1 directory (not build), create mylab.mac:
+
cd ~/g4-labs/lab2/B1
+nano mylab.mac
+
Or use your preferred text editor (gedit, vim, etc.)
+
Add the following commands to your macro file:
+
# My custom run for Lab 2
+# Author: [Your Name]
+
+# Initialize the kernel
+/run/initialize
+
+# Set the particle to a proton
+/gun/particle proton
+
+# Set the proton's kinetic energy to 150 MeV
+/gun/energy 150 MeV
+
+# Print progress every 100 events
+/run/printProgress 100
+
+# Run 1000 events
+/run/beamOn 1000
+
Save and close the file.
+
Run the simulation using your new macro:
+
cd build
+./exampleB1 ../mylab.mac
+
+
+
+
+Custom Macro Execution
+
+
+
+
+
Step 7: Exploring Different Particles and Energies
+
Let’s create additional macro files to explore different simulation scenarios:
+
+
Create a macro for electron simulation (electron_beam.mac):
Compare the output in your terminal from the default run (10 gammas) and your custom run (1000 protons).
+
Look at the final lines of the output for both runs. Find the section that says: The run consists of... The mean energy deposited in the envelope is...
+
How does the mean energy deposited per event differ between:
+
+
6 MeV gamma rays
+
150 MeV protons
+
50 MeV electrons
+
14 MeV neutrons
+
+
Does the result make sense intuitively? Why or why not? Consider:
+
+
Different interaction mechanisms for different particles
+
Energy of the incident particles
+
Range and stopping power
+
+
Why is batch mode more suitable for running 1000 events than interactive mode?
+
What advantages does using macro files provide for scientific reproducibility?
+
+
+
+
Optional: Setting up VS Code IDE
+
For more advanced development work, you may want to use an Integrated Development Environment (IDE) like VS Code:
Error: “Geant4 not found”: Make sure you’ve sourced the geant4.sh script in your current terminal session
+
No GUI appears: Verify that Qt5 was enabled when you built Geant4 (GEANT4_USE_QT=ON)
+
Compilation errors: Ensure all Lab 1 dependencies are installed
+
Macro file not found: Check that you’re running from the build directory and using the correct relative path (../mylab.mac)
+
+
+
+
Verification Checklist
+
Before completing this lab, verify you can:
+
+
+
+
+
+
+
+
+
+
+
Conclusion
+
You have successfully compiled and run your first Geant4 simulation using Example B1. You’ve learned to work in both interactive and batch modes, and you’ve seen how to modify simulation parameters using macro files. These skills form the foundation for all future Geant4 work, where you’ll create increasingly complex simulations for radiation physics and detector design applications.
+
+
+
References
+
+
Geant4 Collaboration. (n.d.). Book for Application Developers
+
Geant4 Collaboration. (n.d.). README for examples/basic/B1. GitLab
+
Geant4 Example B1 Macro Files: run1.mac, exampleB1.in
In this lab, you will bridge the gap between simulation and analysis. You will modify the C++ source code of a Geant4 application to save data and then use Python to analyze and visualize that data. By the end of this session, you will be able to:
+
+
Understand the basic structure of Geant4 user action classes
+
Modify C++ code to implement the G4AnalysisManager for data output
+
Save simulation results to CSV file format for universal compatibility
+
Write Python scripts using the pandas and matplotlib libraries
+
Load simulation data into pandas DataFrames
+
Create and customize histograms to visualize distributions
+
Calculate descriptive statistics (mean, standard deviation, percentiles) from data
+
Filter and manipulate data using pandas operations
+
Compare multiple datasets and create multi-plot figures
+
+
+
+
Introduction
+
A simulation that produces no output is not very useful. To perform scientific analysis, we need to “score” quantities of interest—like the energy deposited in a detector—and save this data to a file. In this lab, you will learn the modern, standard way to do this in Geant4 using a special class called G4AnalysisManager. We will configure it to save the energy deposited in each event into a simple text-based CSV (Comma-Separated Values) file.
+
+
+
+Data Analysis Workflow
+
+
+
Once we have the data, we will switch hats from a simulation user to a data analyst. We will use the powerful Python programming language along with two essential libraries: pandas for data manipulation and matplotlib for plotting. You will write a script to read the CSV file and generate a histogram, which is a graph that shows the frequency distribution of the energy deposition values. This process—simulate, save, analyze, visualize—is a fundamental cycle in computational science.
+
+
+
Prerequisites
+
Before you begin, ensure you have the following:
+
+
Completed Lab 2 (successfully built and ran Example B1)
+
Text editor for C++ code (gedit, nano, or VS Code)
+
Basic understanding of C++ syntax (we’ll guide you through the modifications)
+
Python 3 installed (should be pre-installed on Ubuntu)
+
+
+
+
Pre-Lab Questions
+
+
What is a CSV file? Why is it a good choice for this lab compared to a more complex format like ROOT?
+
What is a histogram and what does it show?
+
Briefly describe the roles of the pandas and matplotlib Python libraries.
+
In Geant4, what is the difference between a “Run”, an “Event”, and a “Step”?
+
What does it mean to “score” a quantity in a simulation?
+
+
+
+
Step 1: Set Up the Lab and Review the Code Structure
+
+
Create a new directory for this lab and copy the B1 example into it, just as you did in Lab 2:
Before we start modifying, let’s understand the key files and the Geant4 action hierarchy:
+
+
+
+Geant4 User Action Classes
+
+
+
In the modern B1 example, all of these classes live inside namespace B1 { ... }. The header and source files do not carry a B1 filename prefix — for example, the run-action class is declared as class RunAction inside namespace B1, in include/RunAction.hh. Every code addition you make in this lab will go inside an existing namespace B1 block in the corresponding .cc file.
+
Familiarize yourself with these key files:
+
+
include/RunAction.hh & src/RunAction.cc:
+
+
Controls actions at the beginning and end of a Run (a collection of events)
+
We will edit the source file to open, write, and close the CSV output
+
+
include/EventAction.hh & src/EventAction.cc:
+
+
Controls actions at the beginning and end of each Event (simulation of one primary particle)
+
We will edit the source file to record one row of CSV data per event
+
+
src/SteppingAction.cc:
+
+
Executed at every single Step a particle takes through the geometry
+
The example already collects per-step energy here. We will not modify this file.
+
+
+
Open src/SteppingAction.cc in a text editor and read through UserSteppingAction:
+
gedit src/SteppingAction.cc &
+
Notice that the example already extracts each step’s energy deposit and adds it to the event total via fEventAction->AddEdep(edepStep). The data lives in memory but never reaches disk — fixing that is what this lab is about.
+
+
+
+
Step 2: Modify C++ Code to Save Data
+
+
What the example already does
+
Before we add anything, look at how energy already flows through the unmodified B1 example. Open these three files in your editor and find the lines below — you do not need to change them:
+
+
include/EventAction.hh already declares the helper and the per-event accumulator:
So the per-step → per-event → per-run pipeline is already wired up. The energy lands in fEdep at the end of each event. Our only job in this lab is to also write that value to a CSV file.
+
+
+
What we’ll change
+
There are exactly two files to edit:
+
+
src/RunAction.cc — set up the CSV file at the start of the run and write/close it at the end.
+
src/EventAction.cc — at the end of each event, write fEdep as one new row of the CSV.
+
+
That’s it. Two #include additions and ~12 new lines of code in total.
+
+
+
+G4AnalysisManager Workflow
+
+
+
+
+
Modification 1: RunAction.cc
+
We use G4AnalysisManager, Geant4’s built-in writer that supports CSV, ROOT, XML, and HDF5. We’ll use CSV.
+
+
Open src/RunAction.cc:
+
gedit src/RunAction.cc &
+
Add the analysis-manager header near the other #include lines at the top:
+
#include "G4AnalysisManager.hh"
+
Inside RunAction::BeginOfRunAction(const G4Run*), add the following block (anywhere inside the function body is fine — putting it at the top is conventional):
+
// --- Lab 3 addition: open CSV and define the per-event ntuple ---
+auto analysisManager = G4AnalysisManager::Instance();
+analysisManager->SetDefaultFileType("csv");
+analysisManager->OpenFile("B1_output");
+
+analysisManager->CreateNtuple("B1","Energy Deposition Data");
+analysisManager->CreateNtupleDColumn("Edep");// column 0
+analysisManager->FinishNtuple();
+
What this does: Tells the analysis manager to write CSV, gives the file a name (B1_output), and declares one ntuple column called Edep. Each row will hold one event’s fEdep.
+
Inside RunAction::EndOfRunAction(const G4Run* run), add the following at the very end of the function (after the existing G4cout print block):
+
// --- Lab 3 addition: flush and close the CSV ---
+auto analysisManager = G4AnalysisManager::Instance();
+analysisManager->Write();
+analysisManager->CloseFile();
+
+
+
+
Modification 2: EventAction.cc
+
+
Open src/EventAction.cc:
+
gedit src/EventAction.cc &
+
Add the analysis-manager header near the top:
+
#include "G4AnalysisManager.hh"
+
Find the existing EventAction::EndOfEventAction function. It currently looks like this — do not retype it, you are only locating it:
+
void EventAction::EndOfEventAction(const G4Event*)
+{
+// accumulate statistics in run action
+ fRunAction->AddEdep(fEdep);
+}
+
Inside that function, after the existing fRunAction->AddEdep(fEdep); line and before the closing }, paste these four lines (and only these four lines):
+
// --- Lab 3 addition: also write fEdep to the CSV ---
+auto analysisManager = G4AnalysisManager::Instance();
+analysisManager->FillNtupleDColumn(0, fEdep);
+analysisManager->AddNtupleRow();
+
After your edit, the complete function should read:
+
void EventAction::EndOfEventAction(const G4Event*)
+{
+// accumulate statistics in run action
+ fRunAction->AddEdep(fEdep);
+
+// --- Lab 3 addition: also write fEdep to the CSV ---
+auto analysisManager = G4AnalysisManager::Instance();
+ analysisManager->FillNtupleDColumn(0, fEdep);
+ analysisManager->AddNtupleRow();
+}
+
What this does: Each event ends with one new row in the CSV containing the total energy deposited in the scoring volume during that event.
+
+
+
+
Summary of Modifications
+
You’ve now modified two files:
+
+
+
+
+
+
+
+
File
+
Change
+
+
+
+
+
src/RunAction.cc
+
Add G4AnalysisManager.hh; open the CSV in BeginOfRunAction, write/close it in EndOfRunAction
+
+
+
src/EventAction.cc
+
Add G4AnalysisManager.hh; in EndOfEventAction, fill one ntuple row per event
+
+
+
+
Note that EventAction.hh, SteppingAction.cc, and the energy-accumulation logic itself were left untouched — the example already had everything we needed up to the point of writing to disk.
+
+
+
+
Step 3: Compile and Run the Modified Application
+
Now let’s compile your modified code and generate some data!
+
+
Navigate to the build directory (create one if it doesn’t exist):
+
mkdir-p build
+cd build
+
Configure the build with CMake:
+
cmake ..
+
Compile the modified application:
+
make-j$(nproc)
+
If you encounter errors, carefully check your code modifications against the instructions. Common issues:
+
+
Missing semicolons
+
Incorrect placement of code
+
Typos in variable names
+
+
+
+
+
+Successful Compilation
+
+
+
+
Create a macro file for data collection. In the B1 directory (not build), create lab3.mac:
+
cd ..
+nano lab3.mac
+
Add these contents:
+
# Lab 3 Data Collection Run
+# Simulate 10,000 events with default 6 MeV gammas
+
+# Force single-threaded mode so we get ONE CSV file, not one per worker.
+# Must come BEFORE /run/initialize.
+/run/numberOfThreads 1
+
+/run/initialize
+/run/printProgress 1000
+/run/beamOn 10000
+
This will run 10,000 events and print progress every 1,000 events. /run/numberOfThreads 1 keeps Geant4 from running multi-threaded, which would otherwise produce one CSV file per worker thread.
+
Run the simulation:
+
cd build
+./exampleB1 ../lab3.mac
+
After completion, list the CSV files:
+
ls-lh*.csv
+
You should see one file named B1_output_nt_B1_t0.csv. Geant4’s CSV writer always appends a _tN worker-thread suffix; with one thread, that’s _t0. If you see multiple _tN.csv files, you forgot /run/numberOfThreads 1 — delete them and re-run.
+
View the first few lines:
+
head-20 B1_output_nt_B1_t0.csv
+
You should see a header row followed by numerical energy values in MeV.
+
+
+
+
+CSV Output Preview
+
+
+
+
+
Step 4: Install Python Analysis Libraries
+
Before we can analyze our data, we need to install the required Python packages.
Why the --break-system-packages flag? Recent versions of Ubuntu (24.04+) follow PEP 668 and mark the system Python environment as “externally managed”, so a plain pip3 install pandas will refuse to run with an error: externally-managed-environment message. The flag tells pip to install into the user site-packages directory anyway, which is fine for this VM-based course where we are not maintaining the host system. (In a production environment you would create a virtual environment with python3 -m venv instead.)
+
These packages provide:
+
+
pandas: Data manipulation and analysis
+
matplotlib: Creating visualizations and plots
+
numpy: Numerical computing (dependency for pandas)
Step 5: Analyze the Data with Python - Basic Analysis
+
Now let’s create a Python script to analyze our simulation data.
+
In your B1/build directory, create analyze_basic.py:
+
nano analyze_basic.py
+
Paste the following content. Important: the code must start at column 0 (no leading whitespace), or Python will reject it with IndentationError: unexpected indent.
Run simulations and rename output files. Remember: each run produces B1_output_nt_B1_t0.csv, so we rename it after each run before it gets overwritten.
+
cd build
+
+# Run gamma simulation (already done, but let's be explicit)
+./exampleB1 ../lab3.mac
+mv B1_output_nt_B1_t0.csv gamma_6MeV.csv
+
+# Run proton simulation
+./exampleB1 ../proton_beam.mac
+mv B1_output_nt_B1_t0.csv proton_150MeV.csv
+
+# Run electron simulation
+./exampleB1 ../electron_beam.mac
+mv B1_output_nt_B1_t0.csv electron_50MeV.csv
+
+
Create the comparison script:
+
nano compare_particles.py
+
Paste the following content (must start at column 0 — no leading whitespace):
Now that you have generated and analyzed the data, answer these questions:
+
+
Part 1: Understanding the Distribution
+
+
Describe the shape of the energy deposition histogram for gamma rays. Is the energy deposition the same for every event?
+
Why is there a distribution of values instead of a single value? Consider:
+
+
Stochastic nature of particle interactions
+
Different interaction mechanisms (photoelectric, Compton, pair production)
+
Variation in secondary particle production
+
+
The peak of the distribution represents the most probable energy deposition. What is this value approximately? How does it compare to the initial energy of the gamma particle (6 MeV)?
+
Explain why the deposited energy is typically lower than the incident energy. (Hint: Does the particle always deposit all of its energy in the detector?)
+
+
+
+
Part 2: Understanding Zero Energy Events
+
+
What percentage of events have zero energy deposition? What does this represent physically?
+
Why would a 6 MeV gamma ray pass through the detector without depositing any energy?
+
+
+
+
Part 3: Comparing Particle Types
+
+
Compare the mean energy deposition for:
+
+
6 MeV gammas
+
150 MeV protons
+
+
50 MeV electrons
+
+
Which particle type deposits the most energy on average? Why?
+
Which particle type has the most consistent energy deposition (smallest standard deviation)? Explain this in terms of the physics of how each particle type interacts with matter.
+
Notice the shape differences in the distributions:
+
+
Gammas: Broad, exponential-like distribution
+
Protons: Sharp, peaked distribution
+
Electrons: Intermediate characteristics
+
+
Explain these shape differences based on the interaction mechanisms of each particle type.
+
+
+
+
Part 4: Pandas Data Analysis
+
+
In the advanced analysis script, we filtered events with data[data['Edep'] > 4.0]. Explain in plain English what this pandas operation does.
+
What is the advantage of using pandas DataFrames versus manipulating the CSV file directly with basic Python?
+
The cumulative distribution plot shows what percentage of events have energy deposition below a certain value. What is the energy value at which 50% of events fall below? What is this called statistically?
+
+
+
+
Part 5: Scientific Workflow
+
+
Why is CSV a good intermediate format for this analysis workflow (Geant4 → CSV → Python)?
+
How would you modify the code to save additional information, such as:
+
+
The number of steps per event
+
Energy deposited in different detector volumes
+
Particle type information for secondary particles
+
+
+
+
+
+
Troubleshooting
+
Common issues and solutions:
+
+
“No module named ‘pandas’”: Install with pip3 install --break-system-packages pandas matplotlib numpy.
+
error: externally-managed-environment from pip: Ubuntu 24.04+ marks the system Python as managed (PEP 668). Re-run with --break-system-packages, as shown in Step 4.
+
Many B1_output_nt_B1_tN.csv files instead of one _t0.csv: Geant4 ran multi-threaded and each worker thread wrote its own CSV. Add /run/numberOfThreads 1 to your macro before/run/initialize, delete the per-thread files (rm B1_output_nt_B1_t*.csv), and re-run.
+
CSV file not found: Make sure you’re running Python from the build directory and that the filename in pd.read_csv(...) matches what ls *.csv shows.
+
Empty histogram: Check that your CSV file has data and proper column headers
+
Compilation errors: Verify all code modifications carefully match the instructions
+
G4AnalysisManager errors: Ensure you’ve included the header file in all necessary source files
+
+
+
+
Verification Checklist
+
Before completing this lab, verify you can:
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Conclusion
+
You have successfully modified a Geant4 application to output data, run simulations with different particle types, and performed sophisticated data analysis using Python’s pandas and matplotlib libraries. You’ve learned:
+
+
How Geant4’s action classes work together to record simulation data
+
The modern approach to data output using G4AnalysisManager
+
Essential pandas operations for data manipulation and analysis
+
Creating publication-quality visualizations with matplotlib
+
The scientific workflow from simulation to analysis to visualization
+
+
These skills are fundamental for computational physics and will be used throughout the remainder of this course and in professional scientific computing work.
Add More Columns: Modify the code to save additional data:
+
+
Track length
+
Number of steps
+
Particle position at end of event
+
+
Different Geometries: Modify the detector geometry in DetectorConstruction.cc and observe how it affects energy deposition
+
Energy Scan: Create a script that runs simulations at multiple energies (e.g., 1, 2, 4, 6, 8, 10 MeV) and plots how mean energy deposition varies with incident energy
+
Interactive Plotting: Modify the Python scripts to use plt.show() for interactive exploration of the data
+
Statistical Tests: Use scipy to perform statistical tests comparing different datasets (t-tests, KS tests, etc.)
Lab 4: Medical Physics Application - A PET Scanner
+
+
+
+
+
+
+
+
Author
+
+
John Murphy and Steven Carlson
+
+
+
+
+
Published
+
+
November 11, 2025
+
+
+
+
+
+
+
+
+
+
+
+
+
Learning Objectives
+
In this final lab, you will apply all the skills you have developed to a more complex and realistic application: a Positron Emission Tomography (PET) scanner. By the end of this session, you will be able to:
+
+
Explain the basic operating principle of a PET scanner
+
Compile and visualize a complex multi-volume detector geometry
+
Use Geant4’s primitive scorer / hits-collection system to record per-crystal energy deposits
+
Generate and interpret an energy spectrum from a simulated radiation detector
+
Identify a “photopeak” and connect it to the underlying physics
+
Reconstruct the spatial arrangement of detector hits and identify coincident events
+
Draw “lines of response” connecting coincident hits, the building block of PET image reconstruction
+
+
+
+
Introduction
+
Positron Emission Tomography (PET) is a medical imaging technique used to visualize metabolic processes in the body, especially in oncology. The patient is given a radiotracer such as Fluorine-18 (F-18), a positron-emitting isotope. When an F-18 nucleus decays, it emits a positron. The positron travels a very short distance through tissue before encountering an electron, and the matter–antimatter pair annihilates, converting their rest mass into two gamma-ray photons.
+
+
+
+Positron Annihilation Diagram
+
+
+
To conserve momentum, the two photons fly off in nearly opposite directions, and each has an energy of 511 keV — the rest mass energy of an electron. A PET scanner is essentially a ring of scintillator detectors designed to detect these photon pairs arriving at almost the same time (in “coincidence”). By recording many such pairs and tracing the line connecting each pair back through the patient, a computer reconstructs an image of where the tracer accumulated.
+
In this lab, we use Geant4 Example B3a, a simplified PET scanner simulation. Our goal is to:
+
+
Record the energy deposited in each scintillator crystal for every event.
+
Plot an energy spectrum and identify the 511 keV photopeak — the signature of a gamma ray fully absorbed by a single crystal.
+
Find events with two photopeak hits in coincidence and draw the lines of response (LOR) that connect them.
+
+
+
+
Prerequisites
+
+
Completed Lab 3 (familiarity with G4AnalysisManager, pandas, and matplotlib)
+
Basic understanding of beta-plus decay and electron-positron annihilation
+
C++ and Python environments configured from previous labs
+
+
+
+
Pre-Lab Questions
+
+
What is positron-electron annihilation, and why does it produce two photons rather than one?
+
What is the energy of each annihilation photon, and why is this specific energy important for PET?
+
In a detector energy spectrum, what does the “photopeak” represent? What is the “Compton continuum”?
+
Why does PET imaging rely on detecting two photons in coincidence rather than one at a time?
+
+
+
+
Step 1: Set Up and Compile the B3a Example
+
The B3 example contains two sub-projects, B3a and B3b. They share geometry and physics, but differ in how event-level statistics are accumulated. We will use B3a.
+
+
Create a lab directory and copy the B3 source tree into it:
Create a build directory and configure with CMake:
+
mkdir build
+cd build
+cmake ..
+
Compile the application:
+
make-j$(nproc)
+
The executable will be named exampleB3a.
+
+
+
+
Step 2: Visualize the PET Scanner Geometry
+
Before running a long simulation, look at the detector to understand what we are simulating.
+
+
Run the application in interactive mode:
+
./exampleB3a
+
The Qt GUI will open. The geometry is more complex than B1: rings of scintillator crystals surround a central cylinder representing a patient phantom.
+
+
+
+B3 PET Scanner Geometry
+
+
+
Fire a few events to see particles produced by F-18 decay:
+
/run/beamOn 5
+
You should see positron tracks (short, in the phantom) followed by pairs of gamma tracks heading to opposite sides of the ring.
+
Rotate the view to confirm the ring geometry, then exit the GUI:
+
exit
+
+
+
+
Step 3: Modify the C++ Code to Record Crystal Hits
+
B3a already uses Geant4’s primitive scorer system. In DetectorConstruction.cc, a G4MultiFunctionalDetector named crystal is attached to each crystal volume, with a G4PSEnergyDeposit primitive scorer named edep. There is also a similar scorer for the patient phantom (patient/dose).
+
What you may not realise is that B3a’s exampleB3a.cc already creates a G4ScoreNtupleWriter and calls SetNtupleMerging(true) on it. This means the framework will automatically write a per-event ntuple for each primitive scorer — including event ID, crystal cell ID, and energy deposit — as soon as we hand it an open output file. We don’t have to write our own EndOfEventAction loop, define our own ntuple, or fill any rows ourselves.
+
In modern Geant4, B3a’s action classes (RunAction, EventAction, ActionInitialization) live inside namespace B3a, and the shared geometry/physics classes live inside namespace B3. The header and source files do not carry a B3a filename prefix. Every code addition you make below will go inside an existing namespace B3a { ... } block.
+
So the entire C++ change in this lab is six new lines in one file: tell the analysis manager to use CSV, open an output file in BeginOfRunAction, and write/close it in EndOfRunAction. That’s it.
+
+
Modification: RunAction.cc
+
+
Open src/RunAction.cc:
+
gedit src/RunAction.cc &
+
Add the analysis-manager header at the top, near the other #include lines:
+
#include "G4AnalysisManager.hh"
+
Inside RunAction::BeginOfRunAction(const G4Run*), add at the top of the function body:
+
// --- Lab 4 addition: route the framework's score-ntuple output to CSV
+auto analysisManager = G4AnalysisManager::Instance();
+analysisManager->SetDefaultFileType("csv");
+analysisManager->OpenFile("B3_output");
+
Inside RunAction::EndOfRunAction(const G4Run*), add at the very end of the function:
+
// --- Lab 4 addition: flush and close the CSV
+auto analysisManager = G4AnalysisManager::Instance();
+analysisManager->Write();
+analysisManager->CloseFile();
+
+
That is the only file you need to edit. EventAction.cc is left untouched — the framework’s G4ScoreNtupleWriter fills the ntuples automatically once a file is open.
+
+
+
+
Step 4: Compile, Run, and Inspect the Output
+
+
Re-compile your modified application:
+
cd build
+make-j$(nproc)
+
If you encounter errors, double-check brace placement, semicolons, and that your additions are inside the correct function.
+
Create a macro file lab4.mac in the B3a directory (one level up from build):
+
cd ..
+nano lab4.mac
+
Add:
+
# Lab 4: PET scanner data run
+
+# Force single-threaded mode so we get ONE CSV file, not one per worker.
+# Must come BEFORE /run/initialize. (Same trick as Lab 3.)
+/run/numberOfThreads 1
+
+/run/initialize
+/run/printProgress 5000
+/run/beamOn 50000
+
Run the simulation in batch mode:
+
cd build
+./exampleB3a ../lab4.mac
+
This will take a few minutes depending on your hardware.
+
Verify the output:
+
ls-lh*.csv
+head B3_output_nt_crystal_edep_t0.csv
+
You will see several CSV files. The framework writes one ntuple per primitive scorer:
+
+
B3_output_nt_crystal_edep_t0.csv — the file we want. Three columns: crystal_edep_eventId, crystal_edep_cell (the crystal copy number), and crystal_edep_score (the energy deposit in MeV).
+
B3_output_nt_patient_dose_t0.csv — patient-phantom dose data. We don’t use it in this lab.
+
The same names without the _t0 suffix are master-thread placeholders containing only headers; you can ignore (or rm) them.
+
+
The first line of B3_output_nt_crystal_edep_t0.csv after the #-prefixed metadata should look like 0,1,0.510999 — event 0, crystal 1, ~511 keV deposited (the photopeak signal we are about to plot).
+
+
+
+
Step 5: Plot the Energy Spectrum
+
Now we switch to Python for analysis. We will plot the distribution of Edep across all hits and look for the 511 keV photopeak.
+
In your B3a/build directory, create analyze_spectrum.py:
+
nano analyze_spectrum.py
+
Paste the following content (must start at column 0 — no leading whitespace):
+
import pandas as pd
+import matplotlib.pyplot as plt
+
+# Geant4's CSV begins with metadata lines starting with '#'. We skip those
+# and supply our own column names matching the C++ ntuple definition.
+data = pd.read_csv('B3_output_nt_crystal_edep_t0.csv',
+ comment='#', header=None,
+ names=['EventID', 'CrystalID', 'Edep'])
+data['Edep_keV'] = data['Edep'] *1000.0
+
+print(f"Total hits recorded: {len(data)}")
+print(f"Unique events with at least one hit: {data['EventID'].nunique()}")
+print(f"Mean Edep: {data['Edep_keV'].mean():.2f} keV")
+print(f"Median Edep: {data['Edep_keV'].median():.2f} keV")
+
+plt.figure(figsize=(11, 6))
+plt.hist(data['Edep_keV'], bins=350, range=(0, 700),
+ color='steelblue', alpha=0.8, edgecolor='black', linewidth=0.3)
+plt.axvline(511, color='red', linestyle='--', linewidth=1.5,
+ label='Expected photopeak: 511 keV')
+
+plt.title('PET Detector Crystal Energy Spectrum')
+plt.xlabel('Energy Deposited per Hit (keV)')
+plt.ylabel('Counts')
+plt.legend()
+plt.grid(True, alpha=0.3)
+plt.tight_layout()
+plt.savefig('pet_spectrum.png', dpi=200)
+print("Spectrum saved to pet_spectrum.png")
+
Run it:
+
python3 analyze_spectrum.py
+eog pet_spectrum.png
+
You should see a sharp peak near 511 keV (the photopeak) sitting on top of a broader, lower-energy distribution (the Compton continuum).
+
+
+
+PET Energy Spectrum
+
+
+
+
+
Step 6: Coincidence Analysis and Lines of Response
+
A single 511 keV hit tells us nothing about where the positron decayed. The PET imaging signal comes from coincidence: two photopeak hits, in the same event, on opposite sides of the detector ring. The line connecting the two hit crystals — the line of response (LOR) — passes through the decay point.
+
We did not record (x, y) directly in the C++ code, but we recorded the crystal copy number, which is enough: B3a’s crystals are arranged in a regular ring, so we can compute each crystal’s position from its copy number. The default geometry in DetectorConstruction.cc lays out:
+
+
32 crystals per ring, evenly spaced around 360°
+
9 axial rings stacked along the scanner’s z-axis
+
Crystals are placed at a fixed inner radius (≈ 32 cm; see the constructor for the exact value)
+
+
Copy numbers run sequentially: copy n belongs to ring n // 32, position n % 32 within that ring.
+
Create analyze_coincidence.py in the build directory:
+
nano analyze_coincidence.py
+
Paste the following content (must start at column 0 — no leading whitespace):
+
import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+
+# --- B3a default geometry constants ---
+# If you have changed the detector, update these to match.
+N_CRYST_PER_RING =32
+N_RINGS =9
+RING_RADIUS_MM =320.0# approximate; verify in DetectorConstruction.cc
+AXIAL_PITCH_MM =30.0# approximate ring-to-ring spacing
+
+
+def crystal_position(crystal_id):
+"""Map a B3a crystal copy number to its approximate (x, y, z) center, in mm."""
+ ring_idx = crystal_id // N_CRYST_PER_RING
+ in_ring = crystal_id % N_CRYST_PER_RING
+ theta = (in_ring +0.5) *2.0* np.pi / N_CRYST_PER_RING
+ x = RING_RADIUS_MM * np.cos(theta)
+ y = RING_RADIUS_MM * np.sin(theta)
+ z = (ring_idx - (N_RINGS -1) /2.0) * AXIAL_PITCH_MM
+return x, y, z
+
+
+# --- Load data and add positions ---
+data = pd.read_csv('B3_output_nt_crystal_edep_t0.csv',
+ comment='#', header=None,
+ names=['EventID', 'CrystalID', 'Edep'])
+data['Edep_keV'] = data['Edep'] *1000.0
+xyz = data['CrystalID'].apply(crystal_position)
+data['X'] = xyz.apply(lambda t: t[0])
+data['Y'] = xyz.apply(lambda t: t[1])
+data['Z'] = xyz.apply(lambda t: t[2])
+
+# --- Hit-position scatter (top view) ---
+plt.figure(figsize=(7, 7))
+plt.scatter(data['X'], data['Y'], s=2, alpha=0.1)
+plt.title('Crystal Hit Positions (Top View)')
+plt.xlabel('X (mm)')
+plt.ylabel('Y (mm)')
+plt.axis('equal')
+plt.grid(True, alpha=0.3)
+plt.tight_layout()
+plt.savefig('hit_positions.png', dpi=200)
+
+# --- Filter for photopeak hits ---
+PHOTOPEAK_LOW =450# keV
+PHOTOPEAK_HIGH =550# keV
+photopeak = data[(data['Edep_keV'] >= PHOTOPEAK_LOW) &
+ (data['Edep_keV'] <= PHOTOPEAK_HIGH)]
+print(f"Hits in photopeak window [{PHOTOPEAK_LOW}, {PHOTOPEAK_HIGH}] keV: "
+f"{len(photopeak)} ({100*len(photopeak)/len(data):.1f}% of all hits)")
+
+# --- Find true coincidences: events with exactly 2 photopeak hits ---
+counts_per_event = photopeak.groupby('EventID').size()
+coincidence_events = counts_per_event[counts_per_event ==2].index
+print(f"Total events processed: {data['EventID'].nunique()}")
+print(f"True coincidences found: {len(coincidence_events)}")
+
+coincident_hits = photopeak[photopeak['EventID'].isin(coincidence_events)]
+
+# --- Draw lines of response ---
+plt.figure(figsize=(8, 8))
+plt.scatter(data['X'], data['Y'], c='lightgray', s=1, alpha=0.1,
+ label='all hits')
+
+plotted =0
+max_lines =200
+for event_id in coincidence_events:
+if plotted >= max_lines:
+break
+ pair = coincident_hits[coincident_hits['EventID'] == event_id]
+ plt.plot(pair['X'].values, pair['Y'].values,
+'r-', alpha=0.3, linewidth=0.6)
+ plotted +=1
+
+plt.title(f'Lines of Response (first {plotted} coincidences)')
+plt.xlabel('X (mm)')
+plt.ylabel('Y (mm)')
+plt.axis('equal')
+plt.grid(True, alpha=0.3)
+plt.tight_layout()
+plt.savefig('lines_of_response.png', dpi=200)
+print("LOR figure saved to lines_of_response.png")
You should see a dense cluster of red lines crossing near the center of the ring — that intersection region marks the location of the simulated source.
+
+
+
+Lines of Response
+
+
+
+
+
Analysis Questions
+
+
Energy Spectrum
+
+
Open pet_spectrum.png. At what energy is the main peak centered? How close is this to the theoretical 511 keV?
+
Why is the photopeak so much taller than the rest of the spectrum?
+
The broad, lower-energy distribution is the Compton continuum — gamma rays that scattered inside a crystal and escaped, depositing only part of their energy. Why is it important for a PET scanner to distinguish photopeak hits from Compton-continuum hits?
+
+
+
+
Coincidence Logic
+
+
We defined a true coincidence as “exactly two hits in the photopeak window in the same event.” What might cause an event to have three photopeak hits? What about one?
+
We used a tight window of 450–550 keV. What would happen to the fraction of accepted events — and to the spatial accuracy of the LORs — if we widened the window to 200–700 keV? (Hint: think about Compton-scattered photons from one crystal landing in another.)
+
+
+
+
Lines of Response
+
+
Look at lines_of_response.png. Do the lines all cross at a single point, or is there a spread? What physical effects limit the spatial resolution of a real PET scanner? (Consider: positron range before annihilation, non-collinearity of the two photons, finite crystal size.)
+
Where in the image is the densest crossing region? What does this tell you about where the simulated source is located?
+
+
+
+
Workflow
+
+
We used the existing primitive scorer rather than writing a custom sensitive detector. What did this save us, and what did we give up (e.g., what information about the hit was not available)?
+
How would you extend the simulation to image two separate point sources in the phantom? What changes in the C++? What changes in the Python?
+
+
+
+
+
Troubleshooting
+
+
Compilation errors mentioning G4AnalysisManager: make sure you added #include "G4AnalysisManager.hh" to RunAction.cc.
+
CSV file not found in Python: confirm the filename in pd.read_csv(...) matches what ls *.csv shows in your build directory.
+
No photopeak in the spectrum: check that you ran enough events (50,000 is a reasonable minimum) and that you are reading the file in keV. Also check the histogram range and binning.
+
Many B3_output_nt_crystal_edep_tN.csv files instead of one _t0.csv: Geant4 ran multi-threaded. Add /run/numberOfThreads 1 to your macro before/run/initialize, delete the per-thread files (rm B3_output_nt_*_t*.csv), and re-run.
+
Empty LOR figure: tighten or relax the photopeak window. With very few simulated events, you may need to lower the lower bound (e.g., 400 keV) to find any coincidences.
+
+
+
+
Verification Checklist
+
Before completing this lab, verify you can:
+
+
+
+
+
+
+
+
+
+
+
+
Conclusion
+
In this lab you simulated a complete PET scanner system. You navigated B3a’s existing primitive-scorer infrastructure, added per-hit data output, and used Python to extract two pieces of physics from the same dataset: the 511 keV photopeak in the energy spectrum, and the lines of response that connect coincident annihilation photons. Together these are the building blocks of PET image reconstruction — combining particle physics (annihilation), detector hardware (scintillators), and software (coincidence logic) to see inside the body.
+
+
+
References
+
+
Geant4 Collaboration. (n.d.). README for examples/basic/B3. GitLab.
Cherry, S. R., Sorenson, J. A., & Phelps, M. E. (2012). Physics in Nuclear Medicine (4th ed.). Elsevier Health Sciences.
+
+
+
+
Optional Extensions
+
+
Energy resolution. Real scintillators broaden the photopeak. Add Gaussian smearing of Edep in your Python analysis (e.g., 10% FWHM at 511 keV) and re-run the spectrum and coincidence steps.
+
3-D LORs. Use the Z column from crystal_position to draw lines of response in 3-D with mpl_toolkits.mplot3d.
+
Two sources. Modify the primary generator to alternate between two source positions, then check whether your LOR plot shows two overlapping clusters of crossings.
+
Custom sensitive detector. Replace the primitive scorer with a custom G4VSensitiveDetector and a CrystalHit class (in namespace B3a) that records position and time directly in C++. Compare the results to the copy-number-based approach used here.