mirror of
https://github.com/jwinterm/geant4-radiation-labs.git
synced 2026-10-17 05:41:25 +02:00
Initial commit: four-lab Geant4 radiation series with Quarto site
This commit is contained in:
jwinterm
+26
@@ -0,0 +1,26 @@
|
||||
# Claude harness
|
||||
.claude/
|
||||
|
||||
# Local working directories with copies of Geant4 source / build artefacts
|
||||
# (do not redistribute Geant4 example sources).
|
||||
lab*files/
|
||||
|
||||
# Third-party PDFs and personal/conference materials
|
||||
literature/
|
||||
reports/
|
||||
|
||||
# Old alternate versions of the labs
|
||||
*_gem.qmd
|
||||
*_gem.pdf
|
||||
|
||||
# Quarto build outputs
|
||||
_site/
|
||||
.quarto/
|
||||
*.html
|
||||
*.pdf
|
||||
|
||||
# OS junk
|
||||
.DS_Store
|
||||
Thumbs.db
|
||||
|
||||
/.quarto/
|
||||
@@ -0,0 +1,32 @@
|
||||
# Geant4 Radiation Labs
|
||||
|
||||
A four-lab series taking you from a fresh Linux install to a working PET-scanner simulation, using [Geant4](https://geant4.web.cern.ch/) for the physics and Python (pandas, matplotlib) for analysis.
|
||||
|
||||
**Read the labs at:** <https://jwinterm.github.io/geant4-radiation-labs/>
|
||||
|
||||
## Labs
|
||||
|
||||
| | Title | Topic |
|
||||
|---|---|---|
|
||||
| 1 | Setting Up Linux and Building Geant4 | Ubuntu VM + Geant4 build with Qt visualization |
|
||||
| 2 | First Simulation: Building and Running B1 | Compile B1, run interactively + via macros, vary particle/energy |
|
||||
| 3 | Acquiring and Analyzing Simulation Data | Patch B1 to write CSV via `G4AnalysisManager`; analyze in Python |
|
||||
| 4 | Medical Physics Application: A PET Scanner | Use B3a PET example; identify 511 keV photopeak; reconstruct lines of response |
|
||||
|
||||
## Building locally
|
||||
|
||||
This site is built with [Quarto](https://quarto.org/). To preview locally:
|
||||
|
||||
```bash
|
||||
quarto preview
|
||||
```
|
||||
|
||||
To regenerate a single lab as PDF:
|
||||
|
||||
```bash
|
||||
quarto render lab1_setup/lab1.qmd --to pdf
|
||||
```
|
||||
|
||||
## Authors
|
||||
|
||||
John Murphy and Steven Carlson
|
||||
+54
@@ -0,0 +1,54 @@
|
||||
project:
|
||||
type: website
|
||||
output-dir: _site
|
||||
render:
|
||||
- "index.qmd"
|
||||
- "lab1_setup/lab1.qmd"
|
||||
- "lab2_exampleB1/lab2.qmd"
|
||||
- "lab3_modifyB1/lab3.qmd"
|
||||
- "lab4_exampleB3/lab4.qmd"
|
||||
|
||||
website:
|
||||
title: "Geant4 Radiation Labs"
|
||||
description: "Hands-on Geant4 simulation labs for medical physics."
|
||||
page-navigation: true
|
||||
navbar:
|
||||
left:
|
||||
- href: index.qmd
|
||||
text: Home
|
||||
- text: "Lab 1: Setup"
|
||||
href: lab1_setup/lab1.qmd
|
||||
- text: "Lab 2: B1"
|
||||
href: lab2_exampleB1/lab2.qmd
|
||||
- text: "Lab 3: Modify B1"
|
||||
href: lab3_modifyB1/lab3.qmd
|
||||
- text: "Lab 4: PET (B3a)"
|
||||
href: lab4_exampleB3/lab4.qmd
|
||||
right:
|
||||
- icon: github
|
||||
href: https://github.com/jwinterm/geant4-radiation-labs
|
||||
aria-label: GitHub repo
|
||||
sidebar:
|
||||
style: "docked"
|
||||
search: true
|
||||
contents:
|
||||
- href: index.qmd
|
||||
text: Introduction
|
||||
- section: "Labs"
|
||||
contents:
|
||||
- lab1_setup/lab1.qmd
|
||||
- lab2_exampleB1/lab2.qmd
|
||||
- lab3_modifyB1/lab3.qmd
|
||||
- lab4_exampleB3/lab4.qmd
|
||||
|
||||
format:
|
||||
html:
|
||||
theme: cosmo
|
||||
toc: true
|
||||
toc-depth: 3
|
||||
toc-location: right
|
||||
code-copy: true
|
||||
code-overflow: wrap
|
||||
code-block-bg: true
|
||||
fig-cap-location: bottom
|
||||
link-external-newwindow: true
|
||||
@@ -0,0 +1,60 @@
|
||||
---
|
||||
title: "Geant4 Radiation Labs"
|
||||
subtitle: "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](https://geant4.web.cern.ch/) 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.
|
||||
|
||||
## The labs
|
||||
|
||||
1. **[Lab 1 — Setting Up Linux and Building Geant4](lab1_setup/lab1.qmd)**
|
||||
Install Ubuntu in a virtual machine and build Geant4 from source with Qt visualisation support.
|
||||
|
||||
2. **[Lab 2 — First Simulation: Building and Running B1](lab2_exampleB1/lab2.qmd)**
|
||||
Compile Geant4's B1 example, run it interactively (Qt GUI) and in batch mode via macros, and modify the particle, energy, and event count.
|
||||
|
||||
3. **[Lab 3 — Acquiring and Analysing Simulation Data](lab3_modifyB1/lab3.qmd)**
|
||||
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.
|
||||
|
||||
4. **[Lab 4 — A PET Scanner](lab4_exampleB3/lab4.qmd)**
|
||||
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
|
||||
|
||||
## Source and PDFs
|
||||
|
||||
The full lab source — `.qmd` files, screenshots, and the Python scripts that generate the conceptual diagrams — lives at <https://github.com/jwinterm/geant4-radiation-labs>.
|
||||
|
||||
PDFs of any lab can be regenerated locally with Quarto:
|
||||
|
||||
```bash
|
||||
quarto render lab1_setup/lab1.qmd --to pdf
|
||||
```
|
||||
|
||||
## Citing
|
||||
|
||||
If this series is useful in your teaching or work, please cite as:
|
||||
|
||||
> Murphy, J. and Carlson, S. (2025). *Geant4 Radiation Labs.* <https://jwinterm.github.io/geant4-radiation-labs/>
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 14 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 122 KiB |
@@ -0,0 +1,198 @@
|
||||
---
|
||||
title: "Lab 1: Setting Up Linux and Building Geant4"
|
||||
author: "John Murphy and Steven Carlson"
|
||||
date: "2025-11-11"
|
||||
---
|
||||
|
||||
# Learning Objectives
|
||||
|
||||
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
|
||||
- VirtualBox or VMware
|
||||
- [Ubuntu 24 ISO Image](https://releases.ubuntu.com/noble/)
|
||||
- Internet connection
|
||||
- At least 50 GB of free disk space and 8 GB of RAM available to the VM
|
||||
|
||||
# Pre-Lab Questions
|
||||
|
||||
1. What is a virtual machine, and why are we using one for this course rather than installing Linux on the host machine directly?
|
||||
2. What is the difference between a *source* distribution of a software package and a *binary* (pre-compiled) distribution?
|
||||
3. What does CMake do? Why do we configure a build with CMake before running `make`?
|
||||
4. What is an environment variable, and why does Geant4 ship a script (`geant4.sh`) for setting them?
|
||||
|
||||
# Step 1: Setting Up the Virtual Machine
|
||||
|
||||
|
||||
|
||||
1. Open VirtualBox or VMware.
|
||||
2. Click "New" to create a virtual machine.
|
||||
3. Name your virtual machine and select Linux type/version (Ubuntu).
|
||||
- Do not use the default VM location (create a new accessible folder instead).
|
||||
4. Copy the ISO file into your chosen accessible folder (e.g., Downloads).
|
||||
5. Select the ISO image file when prompted by VirtualBox.
|
||||
6. Allocate memory (RAM): at least 8 GB.
|
||||
7. Allocate CPU cores (maximum 60% of available cores). Check core count in Windows via `Settings > About Your PC`.
|
||||
8. Create a virtual hard disk (at least 50 GB recommended).
|
||||
9. Start the virtual machine and boot from the ISO.
|
||||
10. 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.
|
||||
11. From the Linux desktop, click "Install Ubuntu" and follow installation prompts to create a permanent installation.
|
||||
12. 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.
|
||||
13. 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)
|
||||
|
||||
1. Open Terminal in Linux (`Ctrl+Alt+T`). Update and install essential packages:
|
||||
|
||||
```bash
|
||||
sudo apt update && sudo apt install -y build-essential \
|
||||
dkms linux-headers-$(uname -r) bzip2 tar
|
||||
```
|
||||
|
||||
2. Insert Guest Additions CD via `Devices > Insert Guest Additions CD image`
|
||||
|
||||
Mount the CD with the following command:
|
||||
|
||||
```bash
|
||||
sudo mount /dev/cdrom /mnt
|
||||
```
|
||||
|
||||
Navigate to that folder with cd (change directory):
|
||||
|
||||
```bash
|
||||
cd /mnt
|
||||
```
|
||||
|
||||
Run:
|
||||
|
||||
```bash
|
||||
sudo ./VBoxLinuxAdditions.run
|
||||
```
|
||||
|
||||
3. Enable fullscreen with `Right Ctrl + F`.
|
||||
4. Enable Shared Clipboard which allows you to paste commands into the terminal:
|
||||
- In VirtualBox, go to `Settings > Advanced` and set Shared Clipboard to "Bidirectional".
|
||||
- Restart the VM.
|
||||
|
||||
# Step 3: Install Required Packages
|
||||
|
||||
Run the following to install the libraries needed to build Geant4 with Qt support:
|
||||
|
||||
```bash
|
||||
sudo apt-get install build-essential libexpat1-dev \
|
||||
libxmu-dev cmake cmake-curses-gui qtbase5-dev qt5-qmake
|
||||
```
|
||||
|
||||
# Step 4: Building Geant4 with Qt and Data
|
||||
|
||||
|
||||
|
||||
1. Download Geant4 source code (tar.gz) from https://geant4.web.cern.ch/download/11.3.2.html.
|
||||
2. Extract and move the folder (`geant4-v11.3.2`) to your home directory.
|
||||
|
||||
3. Create a build folder:
|
||||
|
||||
```bash
|
||||
mkdir -p ~/geant4build
|
||||
cd ~/geant4build
|
||||
```
|
||||
|
||||
4. Configure Geant4 with CMake:
|
||||
|
||||
```bash
|
||||
cmake \
|
||||
-DCMAKE_INSTALL_PREFIX=/usr/local \
|
||||
-DGEANT4_INSTALL_DATA=ON \
|
||||
-DGEANT4_INSTALL_DATADIR=/usr/local/share/Geant4/data \
|
||||
-DGEANT4_USE_QT=ON \
|
||||
~/geant4-v11.3.2
|
||||
```
|
||||
|
||||
5. Compile using all available processors:
|
||||
|
||||
```bash
|
||||
make -j$(nproc)
|
||||
```
|
||||
|
||||
6. Install Geant4:
|
||||
|
||||
```bash
|
||||
sudo make install
|
||||
```
|
||||
|
||||
# Step 5: Post-Installation Setup
|
||||
|
||||
1. Configure library paths:
|
||||
|
||||
```bash
|
||||
echo "/usr/local/lib" | sudo tee \
|
||||
/etc/ld.so.conf.d/geant4.conf
|
||||
sudo ldconfig
|
||||
```
|
||||
|
||||
2. 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:
|
||||
|
||||
```bash
|
||||
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
|
||||
|
||||
1. Verify that the environment is set correctly. In the same terminal where you sourced `geant4.sh`, run:
|
||||
|
||||
```bash
|
||||
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.
|
||||
|
||||
2. Navigate to your installation directory and explore the structure:
|
||||
|
||||
```bash
|
||||
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:
|
||||
|
||||
- [ ] Boot a working Ubuntu virtual machine with shared clipboard and fullscreen support
|
||||
- [ ] Install build dependencies and Qt5 development packages from `apt`
|
||||
- [ ] Run CMake to configure the Geant4 build with Qt and data downloads enabled
|
||||
- [ ] Compile and install Geant4 with `make` and `sudo make install`
|
||||
- [ ] Source `geant4.sh` and confirm `G4INSTALL` and related variables are set
|
||||
- [ ] Locate the Geant4 examples directory under `/usr/local/share`
|
||||
|
||||
## 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.
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 74 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 23 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 98 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 112 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 63 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 182 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 178 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 63 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 139 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 136 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 122 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 154 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 30 KiB |
@@ -0,0 +1,400 @@
|
||||
---
|
||||
title: "Lab 2: First Simulation - Building and Running Geant4 Example B1"
|
||||
author: "John Murphy and Steven Carlson"
|
||||
date: "2025-11-11"
|
||||
---
|
||||
|
||||
## 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.
|
||||
|
||||
|
||||
|
||||
You will learn to run this application in two ways:
|
||||
|
||||
1. **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.
|
||||
2. **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
|
||||
|
||||
1. In the context of a simulation, what is an "event"? What is a "run"?
|
||||
2. What is the purpose of a macro file in Geant4?
|
||||
3. Based on the Example B1 README, what is the default particle and energy being simulated?
|
||||
4. 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.
|
||||
|
||||
1. First, verify where your Geant4 source code is located. If you followed Lab 1, it should be in your home directory:
|
||||
|
||||
```bash
|
||||
ls ~/geant4-v11.3.2/examples/basic/B1
|
||||
```
|
||||
|
||||
2. Create a dedicated directory for lab examples and navigate into it:
|
||||
|
||||
```bash
|
||||
mkdir -p ~/g4-labs/lab2
|
||||
cd ~/g4-labs/lab2
|
||||
```
|
||||
|
||||
3. Copy the source code for Example B1 into your lab directory:
|
||||
|
||||
```bash
|
||||
cp -r ~/geant4-v11.3.2/examples/basic/B1 .
|
||||
```
|
||||
|
||||
4. Navigate into the B1 example directory:
|
||||
|
||||
```bash
|
||||
cd B1
|
||||
```
|
||||
|
||||
5. Explore the contents of the directory:
|
||||
|
||||
```bash
|
||||
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
|
||||
|
||||
|
||||
|
||||
## 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.
|
||||
|
||||
1. First, ensure your Geant4 environment is set up. Open a new terminal and run:
|
||||
|
||||
```bash
|
||||
source /usr/local/bin/geant4.sh
|
||||
```
|
||||
|
||||
Verify the environment variables are set:
|
||||
|
||||
```bash
|
||||
env | grep G4
|
||||
```
|
||||
|
||||
2. Create a build directory for the B1 application inside the B1 folder:
|
||||
|
||||
```bash
|
||||
mkdir build
|
||||
cd build
|
||||
```
|
||||
|
||||
3. Run `cmake` to configure the build. The `..` tells `cmake` to look for the source code in the parent directory:
|
||||
|
||||
```bash
|
||||
cmake ..
|
||||
```
|
||||
|
||||
Because you sourced the `geant4.sh` script, `cmake` automatically knows where to find your Geant4 installation.
|
||||
|
||||
|
||||
|
||||
4. Compile the application using `make`:
|
||||
|
||||
```bash
|
||||
make -j$(nproc)
|
||||
```
|
||||
|
||||
This uses all available CPU cores to speed up compilation.
|
||||
|
||||
5. After compilation finishes, verify the executable was created:
|
||||
|
||||
```bash
|
||||
ls -lh exampleB1
|
||||
```
|
||||
|
||||
You should see the `exampleB1` executable file.
|
||||
|
||||
|
||||
|
||||
## Step 3: Run in Interactive Mode (GUI)
|
||||
|
||||
Now let's run the simulation with the graphical user interface to visualize what's happening.
|
||||
|
||||
1. From the `B1/build` directory, run the executable with no arguments:
|
||||
|
||||
```bash
|
||||
./exampleB1
|
||||
```
|
||||
|
||||
2. 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
|
||||
|
||||
|
||||
|
||||
3. The simulation is now waiting for your commands. Let's fire 10 particles. At the `Idle>` prompt, type:
|
||||
|
||||
```bash
|
||||
/run/beamOn 10
|
||||
```
|
||||
|
||||
4. 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.
|
||||
|
||||
|
||||
|
||||
5. 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
|
||||
|
||||
6. Run more events to see different particle behaviors:
|
||||
|
||||
```bash
|
||||
/run/beamOn 50
|
||||
```
|
||||
|
||||
|
||||
|
||||
7. To exit the application, type `exit` at the prompt:
|
||||
|
||||
```bash
|
||||
exit
|
||||
```
|
||||
|
||||
## Step 4: Understanding the Simulation Output
|
||||
|
||||
1. Look at the terminal output from your interactive run. You should see information including:
|
||||
- Run summary
|
||||
- Mean energy deposited in different volumes
|
||||
- Dose calculations
|
||||
|
||||
|
||||
|
||||
2. 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.
|
||||
|
||||
1. Go back to the main `B1` directory (one level up from `build`):
|
||||
|
||||
```bash
|
||||
cd ..
|
||||
```
|
||||
|
||||
2. View the contents of the provided macro file `run1.mac`:
|
||||
|
||||
```bash
|
||||
cat run1.mac
|
||||
```
|
||||
|
||||
You will see a list of Geant4 commands, including `/run/beamOn 10`.
|
||||
|
||||
|
||||
|
||||
3. Now, run the simulation using the macro file:
|
||||
|
||||
```bash
|
||||
cd build
|
||||
./exampleB1 ../run1.mac
|
||||
```
|
||||
|
||||
4. 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.
|
||||
|
||||
5. Let's examine another macro file that shows more options:
|
||||
|
||||
```bash
|
||||
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.
|
||||
|
||||
1. Create a new macro file for your custom simulation. In the `B1` directory (not `build`), create `mylab.mac`:
|
||||
|
||||
```bash
|
||||
cd ~/g4-labs/lab2/B1
|
||||
nano mylab.mac
|
||||
```
|
||||
|
||||
Or use your preferred text editor (gedit, vim, etc.)
|
||||
|
||||
2. 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
|
||||
```
|
||||
|
||||
3. Save and close the file.
|
||||
|
||||
4. Run the simulation using your new macro:
|
||||
|
||||
```bash
|
||||
cd build
|
||||
./exampleB1 ../mylab.mac
|
||||
```
|
||||
|
||||
|
||||
|
||||
## Step 7: Exploring Different Particles and Energies
|
||||
|
||||
Let's create additional macro files to explore different simulation scenarios:
|
||||
|
||||
1. Create a macro for electron simulation (`electron_beam.mac`):
|
||||
|
||||
```
|
||||
/run/initialize
|
||||
/gun/particle e-
|
||||
/gun/energy 50 MeV
|
||||
/run/beamOn 500
|
||||
```
|
||||
|
||||
2. Create a macro for neutron simulation (`neutron_beam.mac`):
|
||||
|
||||
```
|
||||
/run/initialize
|
||||
/gun/particle neutron
|
||||
/gun/energy 14 MeV
|
||||
/run/beamOn 1000
|
||||
```
|
||||
|
||||
3. Run each simulation and compare the results:
|
||||
|
||||
```bash
|
||||
./exampleB1 ../electron_beam.mac
|
||||
./exampleB1 ../neutron_beam.mac
|
||||
```
|
||||
|
||||
## Analysis Questions
|
||||
|
||||
1. Compare the output in your terminal from the default run (10 gammas) and your custom run (1000 protons).
|
||||
|
||||
2. 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...`
|
||||
|
||||
3. How does the mean energy deposited per event differ between:
|
||||
- 6 MeV gamma rays
|
||||
- 150 MeV protons
|
||||
- 50 MeV electrons
|
||||
- 14 MeV neutrons
|
||||
|
||||
4. 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
|
||||
|
||||
5. Why is batch mode more suitable for running 1000 events than interactive mode?
|
||||
|
||||
6. 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:
|
||||
|
||||
1. Download and install [VS Code](https://code.visualstudio.com/)
|
||||
|
||||
2. Install useful extensions:
|
||||
- C/C++ Extension Pack
|
||||
- CMake Tools
|
||||
- Code Runner
|
||||
|
||||
3. Open the B1 directory in VS Code:
|
||||
|
||||
```bash
|
||||
code ~/g4-labs/lab2/B1
|
||||
```
|
||||
|
||||
4. Explore the source code files to understand how the simulation works:
|
||||
- `exampleB1.cc` - Main program
|
||||
- `src/DetectorConstruction.cc` - Geometry definition
|
||||
- `src/PrimaryGeneratorAction.cc` - Particle gun setup
|
||||
- `src/RunAction.cc` - What happens each run
|
||||
- `src/EventAction.cc` - What happens each event
|
||||
|
||||
|
||||
|
||||
## Troubleshooting
|
||||
|
||||
- **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:
|
||||
|
||||
- [ ] Locate and copy Geant4 example code
|
||||
- [ ] Build a Geant4 application using CMake
|
||||
- [ ] Run simulations in interactive mode with GUI
|
||||
- [ ] Run simulations in batch mode using macro files
|
||||
- [ ] Create custom macro files with different particles and energies
|
||||
- [ ] Interpret the simulation output and energy deposition results
|
||||
|
||||
## 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`
|
||||
- SLAC. (2021). *Geant4 Tutorial: Hands On 1*
|
||||
|
||||
## Additional Resources
|
||||
|
||||
- [Geant4 User's Guide for Application Developers](https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/)
|
||||
- [Geant4 Examples Documentation](https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Examples/examples.html)
|
||||
- [Geant4 Command Reference](https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Control/AllResources/Control/UIcommands/html)
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 149 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 324 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 143 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 285 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 77 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 197 KiB |
@@ -0,0 +1,132 @@
|
||||
"""Generate action_classes.png for Lab 3.
|
||||
|
||||
Shows the nested Geant4 user-action hierarchy (Run -> Event -> Step) as
|
||||
three concentric boxes, each labelled with the corresponding B1 class and
|
||||
the work that class will do after Lab 3's modifications.
|
||||
|
||||
Layout strategy: each layer reserves a "header strip" at its top for its
|
||||
title and 2-3 body lines, and the next inner box starts BELOW that strip
|
||||
so nothing overlaps.
|
||||
"""
|
||||
|
||||
import os
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.patches import FancyBboxPatch
|
||||
|
||||
OUT_DIR = os.path.dirname(os.path.abspath(__file__))
|
||||
OUT_PATH = os.path.join(OUT_DIR, "action_classes.png")
|
||||
|
||||
# Figure coordinates (data units = arbitrary).
|
||||
FIG_W = 11.0
|
||||
FIG_H = 7.4
|
||||
|
||||
# Outer Run box.
|
||||
RUN_BOX = (0.4, 0.6, 10.2, 6.4) # x, y, w, h
|
||||
# Inner Event box. Top is 1.5 units below Run top to leave room for header.
|
||||
EVENT_BOX = (1.05, 1.0, 8.9, 4.55)
|
||||
# Innermost Step box. Top is 1.4 units below Event top.
|
||||
STEP_BOX = (1.7, 1.45, 7.6, 2.7)
|
||||
|
||||
LAYERS = [
|
||||
{
|
||||
"rect": RUN_BOX,
|
||||
"color": "#cfe2ff",
|
||||
"title": "Run — RunAction",
|
||||
"lines": [
|
||||
"BeginOfRunAction(): open CSV file, define ntuple",
|
||||
"EndOfRunAction(): write & close CSV file",
|
||||
],
|
||||
},
|
||||
{
|
||||
"rect": EVENT_BOX,
|
||||
"color": "#fff3cd",
|
||||
"title": "Event — EventAction",
|
||||
"lines": [
|
||||
"BeginOfEventAction(): reset fEdep = 0",
|
||||
"EndOfEventAction(): write fEdep to ntuple",
|
||||
],
|
||||
},
|
||||
{
|
||||
"rect": STEP_BOX,
|
||||
"color": "#d4edda",
|
||||
"title": "Step — SteppingAction",
|
||||
"lines": [
|
||||
"UserSteppingAction(): add",
|
||||
" aStep->GetTotalEnergyDeposit()",
|
||||
"to EventAction's running total",
|
||||
],
|
||||
},
|
||||
]
|
||||
|
||||
|
||||
def main() -> None:
|
||||
fig, ax = plt.subplots(figsize=(FIG_W, FIG_H))
|
||||
ax.set_xlim(0, FIG_W)
|
||||
ax.set_ylim(0, FIG_H)
|
||||
ax.axis("off")
|
||||
|
||||
# Draw boxes back-to-front so inner boxes paint over parent fills.
|
||||
for layer in LAYERS:
|
||||
x, y, w, h = layer["rect"]
|
||||
patch = FancyBboxPatch(
|
||||
(x, y),
|
||||
w,
|
||||
h,
|
||||
boxstyle="round,pad=0.05",
|
||||
facecolor=layer["color"],
|
||||
edgecolor="black",
|
||||
linewidth=1.6,
|
||||
)
|
||||
ax.add_patch(patch)
|
||||
|
||||
# Place title + body in each layer's header strip (above the inner box).
|
||||
for layer in LAYERS:
|
||||
x, y, w, h = layer["rect"]
|
||||
title_y = y + h - 0.45
|
||||
ax.text(
|
||||
x + 0.3,
|
||||
title_y,
|
||||
layer["title"],
|
||||
ha="left",
|
||||
va="center",
|
||||
fontsize=12,
|
||||
fontweight="bold",
|
||||
)
|
||||
for i, line in enumerate(layer["lines"]):
|
||||
ax.text(
|
||||
x + 0.45,
|
||||
title_y - 0.4 - i * 0.32,
|
||||
line,
|
||||
ha="left",
|
||||
va="center",
|
||||
fontsize=10,
|
||||
family="monospace",
|
||||
)
|
||||
|
||||
ax.set_title(
|
||||
"Geant4 User-Action Hierarchy in Example B1\n"
|
||||
"(all classes are in namespace B1)",
|
||||
fontsize=13,
|
||||
fontweight="bold",
|
||||
pad=10,
|
||||
)
|
||||
|
||||
ax.text(
|
||||
FIG_W / 2,
|
||||
0.25,
|
||||
"Each Run contains many Events; each Event contains many Steps.",
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=10,
|
||||
style="italic",
|
||||
color="#444",
|
||||
)
|
||||
|
||||
plt.tight_layout()
|
||||
fig.savefig(OUT_PATH, dpi=200, bbox_inches="tight")
|
||||
print(f"Wrote {OUT_PATH}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,208 @@
|
||||
"""Generate analysis_manager_flow.png for Lab 3.
|
||||
|
||||
Shows how a single event's energy flows from per-step calls in
|
||||
SteppingAction, through per-event accumulation in EventAction, into
|
||||
G4AnalysisManager, and finally out to a CSV file. The diagram also
|
||||
calls out the per-run setup/teardown in RunAction (open/write/close).
|
||||
"""
|
||||
|
||||
import os
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.patches import FancyArrowPatch, FancyBboxPatch
|
||||
|
||||
OUT_DIR = os.path.dirname(os.path.abspath(__file__))
|
||||
OUT_PATH = os.path.join(OUT_DIR, "analysis_manager_flow.png")
|
||||
|
||||
# Four boxes in the main left-to-right flow.
|
||||
MAIN_BOXES = [
|
||||
{
|
||||
"title": "SteppingAction",
|
||||
"scope": "(per step)",
|
||||
"code": [
|
||||
"edepStep = aStep->",
|
||||
" GetTotalEnergyDeposit();",
|
||||
"fEventAction->",
|
||||
" AddEdep(edepStep);",
|
||||
],
|
||||
"color": "#cfe2ff",
|
||||
},
|
||||
{
|
||||
"title": "EventAction",
|
||||
"scope": "(per event)",
|
||||
"code": [
|
||||
"fEdep += edepStep;",
|
||||
"",
|
||||
"// EndOfEventAction:",
|
||||
"FillNtupleDColumn(0, fEdep);",
|
||||
"AddNtupleRow();",
|
||||
],
|
||||
"color": "#fff3cd",
|
||||
},
|
||||
{
|
||||
"title": "G4AnalysisManager",
|
||||
"scope": "(buffered)",
|
||||
"code": [
|
||||
"Buffers each row",
|
||||
"until run ends",
|
||||
],
|
||||
"color": "#e2d4f0",
|
||||
},
|
||||
{
|
||||
"title": "CSV File",
|
||||
"scope": "(on disk)",
|
||||
"code": [
|
||||
"B1_output_nt_B1.csv",
|
||||
"",
|
||||
"Edep",
|
||||
"0.043",
|
||||
"0.512",
|
||||
"...",
|
||||
],
|
||||
"color": "#d4edda",
|
||||
},
|
||||
]
|
||||
|
||||
|
||||
def main() -> None:
|
||||
fig, ax = plt.subplots(figsize=(13, 5.4))
|
||||
ax.set_xlim(0, 13)
|
||||
ax.set_ylim(0, 5.4)
|
||||
ax.axis("off")
|
||||
|
||||
box_w, box_h = 2.7, 3.4
|
||||
gap = 0.45
|
||||
total_w = len(MAIN_BOXES) * box_w + (len(MAIN_BOXES) - 1) * gap
|
||||
x0 = (13 - total_w) / 2
|
||||
y0 = 1.1
|
||||
|
||||
centers = []
|
||||
for i, b in enumerate(MAIN_BOXES):
|
||||
x = x0 + i * (box_w + gap)
|
||||
centers.append((x + box_w / 2, y0 + box_h / 2, x, x + box_w))
|
||||
patch = FancyBboxPatch(
|
||||
(x, y0),
|
||||
box_w,
|
||||
box_h,
|
||||
boxstyle="round,pad=0.08",
|
||||
facecolor=b["color"],
|
||||
edgecolor="black",
|
||||
linewidth=1.5,
|
||||
)
|
||||
ax.add_patch(patch)
|
||||
ax.text(
|
||||
x + box_w / 2,
|
||||
y0 + box_h - 0.45,
|
||||
b["title"],
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=11.5,
|
||||
fontweight="bold",
|
||||
)
|
||||
ax.text(
|
||||
x + box_w / 2,
|
||||
y0 + box_h - 0.85,
|
||||
b["scope"],
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=9,
|
||||
style="italic",
|
||||
color="#444",
|
||||
)
|
||||
for j, line in enumerate(b["code"]):
|
||||
ax.text(
|
||||
x + 0.18,
|
||||
y0 + box_h - 1.30 - j * 0.32,
|
||||
line,
|
||||
ha="left",
|
||||
va="center",
|
||||
fontsize=9,
|
||||
family="monospace",
|
||||
)
|
||||
|
||||
# Arrows between consecutive main boxes.
|
||||
arrow_y = y0 + box_h / 2
|
||||
for (_, _, _, right_a), (_, _, left_b, _) in zip(centers, centers[1:]):
|
||||
arrow = FancyArrowPatch(
|
||||
(right_a, arrow_y),
|
||||
(left_b, arrow_y),
|
||||
arrowstyle="-|>",
|
||||
mutation_scale=20,
|
||||
linewidth=1.6,
|
||||
color="black",
|
||||
)
|
||||
ax.add_patch(arrow)
|
||||
|
||||
# RunAction callout: arrow from below pointing up at the AnalysisManager
|
||||
# box, with an explanatory label.
|
||||
am_cx = centers[2][0]
|
||||
am_left = centers[2][2]
|
||||
am_right = centers[2][3]
|
||||
callout_y_text = 0.38
|
||||
ax.text(
|
||||
am_cx,
|
||||
callout_y_text,
|
||||
"RunAction (per run)",
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=10.5,
|
||||
fontweight="bold",
|
||||
)
|
||||
ax.text(
|
||||
am_cx,
|
||||
callout_y_text - 0.32,
|
||||
"BeginOfRunAction: OpenFile + CreateNtuple",
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=9,
|
||||
family="monospace",
|
||||
)
|
||||
callout_arrow = FancyArrowPatch(
|
||||
(am_left + 0.2, callout_y_text + 0.18),
|
||||
(am_left + 0.2, y0 - 0.05),
|
||||
arrowstyle="-|>",
|
||||
mutation_scale=14,
|
||||
linewidth=1.2,
|
||||
color="#444",
|
||||
linestyle="dashed",
|
||||
)
|
||||
ax.add_patch(callout_arrow)
|
||||
|
||||
# EndOfRunAction arrow from the AM box across to the CSV box.
|
||||
csv_left = centers[3][2]
|
||||
end_arrow = FancyArrowPatch(
|
||||
(am_right - 0.2, y0 - 0.05),
|
||||
(csv_left + 0.2, y0 - 0.05),
|
||||
arrowstyle="-|>",
|
||||
mutation_scale=14,
|
||||
linewidth=1.2,
|
||||
color="#444",
|
||||
linestyle="dashed",
|
||||
connectionstyle="arc3,rad=-0.25",
|
||||
)
|
||||
ax.add_patch(end_arrow)
|
||||
ax.text(
|
||||
(am_right + csv_left) / 2,
|
||||
y0 - 0.55,
|
||||
"EndOfRunAction: Write() + CloseFile()",
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=9,
|
||||
family="monospace",
|
||||
color="#444",
|
||||
)
|
||||
|
||||
ax.set_title(
|
||||
"Energy Data Flow: Step → Event → AnalysisManager → CSV",
|
||||
fontsize=13,
|
||||
fontweight="bold",
|
||||
pad=10,
|
||||
)
|
||||
|
||||
plt.tight_layout()
|
||||
fig.savefig(OUT_PATH, dpi=200, bbox_inches="tight")
|
||||
print(f"Wrote {OUT_PATH}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,109 @@
|
||||
"""Generate workflow_diagram.png for Lab 3.
|
||||
|
||||
Renders the Simulate -> Save -> Analyze -> Visualize cycle as four labelled
|
||||
boxes connected by arrows. Re-run this script after editing to regenerate
|
||||
the PNG.
|
||||
"""
|
||||
|
||||
import os
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.patches import FancyArrowPatch, FancyBboxPatch
|
||||
|
||||
OUT_DIR = os.path.dirname(os.path.abspath(__file__))
|
||||
OUT_PATH = os.path.join(OUT_DIR, "workflow_diagram.png")
|
||||
|
||||
BOXES = [
|
||||
{
|
||||
"title": "Geant4 Simulation",
|
||||
"subtitle": "C++ application\n(modified B1)",
|
||||
"color": "#cfe2ff",
|
||||
},
|
||||
{
|
||||
"title": "CSV File",
|
||||
"subtitle": "B1_output.csv\none row per event",
|
||||
"color": "#fff3cd",
|
||||
},
|
||||
{
|
||||
"title": "Python + pandas",
|
||||
"subtitle": "Load DataFrame\nfilter & summarize",
|
||||
"color": "#d4edda",
|
||||
},
|
||||
{
|
||||
"title": "Plots & Statistics",
|
||||
"subtitle": "matplotlib\nhistograms, CDFs",
|
||||
"color": "#f8d7da",
|
||||
},
|
||||
]
|
||||
|
||||
|
||||
def main() -> None:
|
||||
fig, ax = plt.subplots(figsize=(12, 3.6))
|
||||
ax.set_xlim(0, 12)
|
||||
ax.set_ylim(0, 3.5)
|
||||
ax.axis("off")
|
||||
|
||||
box_w, box_h = 2.5, 1.7
|
||||
gap = 0.4
|
||||
total_w = len(BOXES) * box_w + (len(BOXES) - 1) * gap
|
||||
x0 = (12 - total_w) / 2
|
||||
y0 = 0.85
|
||||
|
||||
centers = []
|
||||
for i, b in enumerate(BOXES):
|
||||
x = x0 + i * (box_w + gap)
|
||||
centers.append((x + box_w / 2, y0 + box_h / 2))
|
||||
patch = FancyBboxPatch(
|
||||
(x, y0),
|
||||
box_w,
|
||||
box_h,
|
||||
boxstyle="round,pad=0.08",
|
||||
facecolor=b["color"],
|
||||
edgecolor="black",
|
||||
linewidth=1.5,
|
||||
)
|
||||
ax.add_patch(patch)
|
||||
ax.text(
|
||||
x + box_w / 2,
|
||||
y0 + box_h * 0.68,
|
||||
b["title"],
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=12,
|
||||
fontweight="bold",
|
||||
)
|
||||
ax.text(
|
||||
x + box_w / 2,
|
||||
y0 + box_h * 0.28,
|
||||
b["subtitle"],
|
||||
ha="center",
|
||||
va="center",
|
||||
fontsize=10,
|
||||
)
|
||||
|
||||
# Arrows between consecutive boxes.
|
||||
for (cx1, cy), (cx2, _) in zip(centers, centers[1:]):
|
||||
arrow = FancyArrowPatch(
|
||||
(cx1 + box_w / 2, cy),
|
||||
(cx2 - box_w / 2, cy),
|
||||
arrowstyle="-|>",
|
||||
mutation_scale=20,
|
||||
linewidth=1.6,
|
||||
color="black",
|
||||
)
|
||||
ax.add_patch(arrow)
|
||||
|
||||
ax.set_title(
|
||||
"Lab 3 Workflow: Simulate → Save → Analyze → Visualize",
|
||||
fontsize=13,
|
||||
fontweight="bold",
|
||||
pad=12,
|
||||
)
|
||||
|
||||
plt.tight_layout()
|
||||
fig.savefig(OUT_PATH, dpi=200, bbox_inches="tight")
|
||||
print(f"Wrote {OUT_PATH}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 313 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 77 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 84 KiB |
@@ -0,0 +1,840 @@
|
||||
---
|
||||
title: "Lab 3: Acquiring and Analyzing Simulation Data"
|
||||
author: "John Murphy and Steven Carlson"
|
||||
date: "2025-11-11"
|
||||
---
|
||||
|
||||
## Learning Objectives
|
||||
|
||||
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.
|
||||
|
||||
|
||||
|
||||
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
|
||||
|
||||
1. What is a CSV file? Why is it a good choice for this lab compared to a more complex format like ROOT?
|
||||
2. What is a histogram and what does it show?
|
||||
3. Briefly describe the roles of the `pandas` and `matplotlib` Python libraries.
|
||||
4. In Geant4, what is the difference between a "Run", an "Event", and a "Step"?
|
||||
5. What does it mean to "score" a quantity in a simulation?
|
||||
|
||||
## Step 1: Set Up the Lab and Review the Code Structure
|
||||
|
||||
1. Create a new directory for this lab and copy the B1 example into it, just as you did in Lab 2:
|
||||
|
||||
```bash
|
||||
mkdir -p ~/g4-labs/lab3
|
||||
cd ~/g4-labs/lab3
|
||||
cp -r ~/geant4-v11.3.2/examples/basic/B1 .
|
||||
cd B1
|
||||
```
|
||||
|
||||
2. Before we start modifying, let's understand the key files and the Geant4 action hierarchy:
|
||||
|
||||
|
||||
|
||||
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.
|
||||
|
||||
3. Open `src/SteppingAction.cc` in a text editor and read through `UserSteppingAction`:
|
||||
|
||||
```bash
|
||||
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:
|
||||
|
||||
```cpp
|
||||
void AddEdep(G4double edep) { fEdep += edep; }
|
||||
G4double fEdep = 0.;
|
||||
```
|
||||
|
||||
- `src/EventAction.cc` already resets `fEdep` at the start of each event and forwards it to the run action at the end:
|
||||
|
||||
```cpp
|
||||
void EventAction::BeginOfEventAction(const G4Event*) { fEdep = 0.; }
|
||||
void EventAction::EndOfEventAction(const G4Event*) { fRunAction->AddEdep(fEdep); }
|
||||
```
|
||||
|
||||
- `src/SteppingAction.cc` already accumulates each step's deposit into the event total:
|
||||
|
||||
```cpp
|
||||
G4double edepStep = step->GetTotalEnergyDeposit();
|
||||
fEventAction->AddEdep(edepStep);
|
||||
```
|
||||
|
||||
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:
|
||||
|
||||
1. **`src/RunAction.cc`** — set up the CSV file at the start of the run and write/close it at the end.
|
||||
2. **`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.
|
||||
|
||||
|
||||
|
||||
### Modification 1: RunAction.cc
|
||||
|
||||
We use `G4AnalysisManager`, Geant4's built-in writer that supports CSV, ROOT, XML, and HDF5. We'll use CSV.
|
||||
|
||||
1. Open `src/RunAction.cc`:
|
||||
|
||||
```bash
|
||||
gedit src/RunAction.cc &
|
||||
```
|
||||
|
||||
2. Add the analysis-manager header near the other `#include` lines at the top:
|
||||
|
||||
```cpp
|
||||
#include "G4AnalysisManager.hh"
|
||||
```
|
||||
|
||||
3. Inside `RunAction::BeginOfRunAction(const G4Run*)`, add the following block (anywhere inside the function body is fine — putting it at the top is conventional):
|
||||
|
||||
```cpp
|
||||
// --- 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`.
|
||||
|
||||
4. Inside `RunAction::EndOfRunAction(const G4Run* run)`, add the following at the **very end** of the function (after the existing `G4cout` print block):
|
||||
|
||||
```cpp
|
||||
// --- Lab 3 addition: flush and close the CSV ---
|
||||
auto analysisManager = G4AnalysisManager::Instance();
|
||||
analysisManager->Write();
|
||||
analysisManager->CloseFile();
|
||||
```
|
||||
|
||||
### Modification 2: EventAction.cc
|
||||
|
||||
1. Open `src/EventAction.cc`:
|
||||
|
||||
```bash
|
||||
gedit src/EventAction.cc &
|
||||
```
|
||||
|
||||
2. Add the analysis-manager header near the top:
|
||||
|
||||
```cpp
|
||||
#include "G4AnalysisManager.hh"
|
||||
```
|
||||
|
||||
3. Find the existing `EventAction::EndOfEventAction` function. It currently looks like this — **do not retype it**, you are only locating it:
|
||||
|
||||
```cpp
|
||||
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):
|
||||
|
||||
```cpp
|
||||
// --- 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:
|
||||
|
||||
```cpp
|
||||
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!
|
||||
|
||||
1. Navigate to the build directory (create one if it doesn't exist):
|
||||
|
||||
```bash
|
||||
mkdir -p build
|
||||
cd build
|
||||
```
|
||||
|
||||
2. Configure the build with CMake:
|
||||
|
||||
```bash
|
||||
cmake ..
|
||||
```
|
||||
|
||||
3. Compile the modified application:
|
||||
|
||||
```bash
|
||||
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
|
||||
|
||||
|
||||
|
||||
4. Create a macro file for data collection. In the `B1` directory (not `build`), create `lab3.mac`:
|
||||
|
||||
```bash
|
||||
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.
|
||||
|
||||
5. Run the simulation:
|
||||
|
||||
```bash
|
||||
cd build
|
||||
./exampleB1 ../lab3.mac
|
||||
```
|
||||
|
||||
6. After completion, list the CSV files:
|
||||
|
||||
```bash
|
||||
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.
|
||||
|
||||
7. View the first few lines:
|
||||
|
||||
```bash
|
||||
head -20 B1_output_nt_B1_t0.csv
|
||||
```
|
||||
|
||||
You should see a header row followed by numerical energy values in MeV.
|
||||
|
||||
|
||||
|
||||
## Step 4: Install Python Analysis Libraries
|
||||
|
||||
Before we can analyze our data, we need to install the required Python packages.
|
||||
|
||||
1. First, check if pip is installed:
|
||||
|
||||
```bash
|
||||
pip3 --version
|
||||
```
|
||||
|
||||
2. If needed, install pip:
|
||||
|
||||
```bash
|
||||
sudo apt install python3-pip
|
||||
```
|
||||
|
||||
3. Install pandas and matplotlib:
|
||||
|
||||
```bash
|
||||
pip3 install --break-system-packages pandas matplotlib numpy
|
||||
```
|
||||
|
||||
**Why the `--break-system-packages` flag?** Recent versions of Ubuntu (24.04+) follow [PEP 668](https://peps.python.org/pep-0668/) 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)
|
||||
|
||||
|
||||
|
||||
4. Verify the installation:
|
||||
|
||||
```bash
|
||||
python3 -c "import pandas; import matplotlib; print('Success!')"
|
||||
```
|
||||
|
||||
## 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`:
|
||||
|
||||
```bash
|
||||
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`.
|
||||
|
||||
```python
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# --- Load the Data ---
|
||||
# Geant4's CSV output begins with metadata lines starting with '#'
|
||||
# (e.g. "#class tools::wcsv::ntuple", "#column double Edep").
|
||||
# We tell pandas to skip those and supply the column name ourselves.
|
||||
print("=" * 50)
|
||||
print("GEANT4 SIMULATION DATA ANALYSIS")
|
||||
print("=" * 50)
|
||||
print("\nReading data from B1_output_nt_B1_t0.csv...")
|
||||
|
||||
data = pd.read_csv('B1_output_nt_B1_t0.csv',
|
||||
comment='#', header=None, names=['Edep'])
|
||||
|
||||
print(f"Data loaded: {len(data)} events, shape {data.shape}")
|
||||
print("\nFirst 10 energy deposition values (MeV):")
|
||||
print(data.head(10))
|
||||
|
||||
# --- Basic Statistics ---
|
||||
print("\n" + "=" * 50)
|
||||
print("BASIC STATISTICS")
|
||||
print("=" * 50)
|
||||
|
||||
mean_edep = data['Edep'].mean()
|
||||
std_edep = data['Edep'].std()
|
||||
median_edep = data['Edep'].median()
|
||||
|
||||
print(f"Mean Energy Deposition: {mean_edep:.4f} MeV")
|
||||
print(f"Standard Deviation: {std_edep:.4f} MeV")
|
||||
print(f"Median: {median_edep:.4f} MeV")
|
||||
print(f"Minimum: {data['Edep'].min():.4f} MeV")
|
||||
print(f"Maximum: {data['Edep'].max():.4f} MeV")
|
||||
|
||||
percentiles = data['Edep'].quantile([0.25, 0.5, 0.75, 0.90, 0.95])
|
||||
print(f"\nPercentiles:")
|
||||
print(f" 25th percentile (Q1): {percentiles[0.25]:.4f} MeV")
|
||||
print(f" 50th percentile (median): {percentiles[0.50]:.4f} MeV")
|
||||
print(f" 75th percentile (Q3): {percentiles[0.75]:.4f} MeV")
|
||||
print(f" 90th percentile: {percentiles[0.90]:.4f} MeV")
|
||||
print(f" 95th percentile: {percentiles[0.95]:.4f} MeV")
|
||||
|
||||
# --- Histogram ---
|
||||
plt.figure(figsize=(12, 8))
|
||||
plt.hist(data['Edep'], bins=100, alpha=0.75,
|
||||
color='steelblue', edgecolor='black')
|
||||
plt.axvline(mean_edep, color='red', linestyle='--',
|
||||
linewidth=2, label=f'Mean: {mean_edep:.3f} MeV')
|
||||
plt.axvline(median_edep, color='green', linestyle='--',
|
||||
linewidth=2, label=f'Median: {median_edep:.3f} MeV')
|
||||
|
||||
plt.title(f'Energy Deposition Distribution for 6 MeV Gammas\n({len(data)} events)',
|
||||
fontsize=14, fontweight='bold')
|
||||
plt.xlabel('Energy Deposited (MeV)', fontsize=12)
|
||||
plt.ylabel('Frequency (Number of Events)', fontsize=12)
|
||||
plt.grid(True, alpha=0.3)
|
||||
plt.legend(fontsize=10)
|
||||
|
||||
stats_text = f'mu = {mean_edep:.3f} MeV\nsigma = {std_edep:.3f} MeV\nEvents = {len(data)}'
|
||||
plt.text(0.98, 0.97, stats_text, transform=plt.gca().transAxes,
|
||||
fontsize=11, verticalalignment='top', horizontalalignment='right',
|
||||
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig('edep_histogram.png', dpi=300)
|
||||
print("\nHistogram saved to edep_histogram.png")
|
||||
```
|
||||
|
||||
Save and close, then run:
|
||||
|
||||
```bash
|
||||
python3 analyze_basic.py
|
||||
```
|
||||
|
||||
View the generated histogram:
|
||||
|
||||
```bash
|
||||
eog edep_histogram.png
|
||||
```
|
||||
|
||||
(Or use your preferred image viewer.)
|
||||
|
||||
|
||||
|
||||
## Step 6: Advanced Analysis - Data Filtering and Comparison
|
||||
|
||||
Let's explore pandas' powerful data manipulation capabilities by filtering the data and creating comparative visualizations.
|
||||
|
||||
Create `analyze_advanced.py`:
|
||||
|
||||
```bash
|
||||
nano analyze_advanced.py
|
||||
```
|
||||
|
||||
Paste the following content (must start at column 0 — no leading whitespace):
|
||||
|
||||
```python
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
print("Loading simulation data...")
|
||||
data = pd.read_csv('B1_output_nt_B1_t0.csv',
|
||||
comment='#', header=None, names=['Edep'])
|
||||
|
||||
# --- Data Filtering Examples ---
|
||||
print("\n" + "=" * 50)
|
||||
print("DATA FILTERING EXAMPLES")
|
||||
print("=" * 50)
|
||||
|
||||
high_energy = data[data['Edep'] > 4.0]
|
||||
print(f"\nEvents with Edep > 4.0 MeV: {len(high_energy)} "
|
||||
f"({len(high_energy)/len(data)*100:.1f}%)")
|
||||
|
||||
low_energy = data[data['Edep'] < 1.0]
|
||||
print(f"Events with Edep < 1.0 MeV: {len(low_energy)} "
|
||||
f"({len(low_energy)/len(data)*100:.1f}%)")
|
||||
|
||||
zero_energy = data[data['Edep'] == 0.0]
|
||||
print(f"Events with Edep = 0.0 MeV: {len(zero_energy)} "
|
||||
f"({len(zero_energy)/len(data)*100:.1f}%)")
|
||||
print(" -> These are events where the particle passed through "
|
||||
"without interacting")
|
||||
|
||||
mid_range = data[(data['Edep'] >= 2.0) & (data['Edep'] <= 4.0)]
|
||||
print(f"Events with 2.0 <= Edep <= 4.0 MeV: {len(mid_range)} "
|
||||
f"({len(mid_range)/len(data)*100:.1f}%)")
|
||||
|
||||
# --- Multi-Panel Figure ---
|
||||
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
|
||||
fig.suptitle('Advanced Energy Deposition Analysis',
|
||||
fontsize=16, fontweight='bold')
|
||||
|
||||
# Panel 1: Full distribution
|
||||
ax1 = axes[0, 0]
|
||||
ax1.hist(data['Edep'], bins=100, color='steelblue',
|
||||
alpha=0.7, edgecolor='black')
|
||||
ax1.set_xlabel('Energy Deposited (MeV)')
|
||||
ax1.set_ylabel('Frequency')
|
||||
ax1.set_title(f'Full Distribution (n={len(data)})')
|
||||
ax1.grid(True, alpha=0.3)
|
||||
|
||||
# Panel 2: Excluding zero-energy events
|
||||
ax2 = axes[0, 1]
|
||||
non_zero = data[data['Edep'] > 0]
|
||||
ax2.hist(non_zero['Edep'], bins=80, color='coral',
|
||||
alpha=0.7, edgecolor='black')
|
||||
ax2.set_xlabel('Energy Deposited (MeV)')
|
||||
ax2.set_ylabel('Frequency')
|
||||
ax2.set_title(f'Non-Zero Events Only (n={len(non_zero)})')
|
||||
ax2.grid(True, alpha=0.3)
|
||||
|
||||
# Panel 3: Cumulative distribution
|
||||
ax3 = axes[1, 0]
|
||||
sorted_data = np.sort(data['Edep'])
|
||||
cumulative = np.arange(1, len(sorted_data) + 1) / len(sorted_data) * 100
|
||||
ax3.plot(sorted_data, cumulative, linewidth=2, color='darkgreen')
|
||||
ax3.set_xlabel('Energy Deposited (MeV)')
|
||||
ax3.set_ylabel('Cumulative Percentage (%)')
|
||||
ax3.set_title('Cumulative Distribution Function')
|
||||
ax3.grid(True, alpha=0.3)
|
||||
ax3.axhline(50, color='red', linestyle='--', alpha=0.5, label='50%')
|
||||
ax3.axhline(95, color='orange', linestyle='--', alpha=0.5, label='95%')
|
||||
ax3.legend()
|
||||
|
||||
# Panel 4: Box plot
|
||||
ax4 = axes[1, 1]
|
||||
box_data = [data['Edep'], non_zero['Edep']]
|
||||
bp = ax4.boxplot(box_data, labels=['All Events', 'Non-Zero Only'],
|
||||
patch_artist=True, showmeans=True)
|
||||
for patch, color in zip(bp['boxes'], ['lightblue', 'lightcoral']):
|
||||
patch.set_facecolor(color)
|
||||
ax4.set_ylabel('Energy Deposited (MeV)')
|
||||
ax4.set_title('Box Plot Comparison')
|
||||
ax4.grid(True, alpha=0.3, axis='y')
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig('edep_advanced_analysis.png', dpi=300)
|
||||
print("\nAdvanced analysis figure saved to edep_advanced_analysis.png")
|
||||
|
||||
# --- Statistical Comparison ---
|
||||
print("\n" + "=" * 50)
|
||||
print("STATISTICAL COMPARISON")
|
||||
print("=" * 50)
|
||||
print("\nAll Events:")
|
||||
print(f" Mean: {data['Edep'].mean():.4f} MeV")
|
||||
print(f" Median: {data['Edep'].median():.4f} MeV")
|
||||
print(f" Std: {data['Edep'].std():.4f} MeV")
|
||||
print("\nNon-Zero Events Only:")
|
||||
print(f" Mean: {non_zero['Edep'].mean():.4f} MeV")
|
||||
print(f" Median: {non_zero['Edep'].median():.4f} MeV")
|
||||
print(f" Std: {non_zero['Edep'].std():.4f} MeV")
|
||||
|
||||
print("\n" + "=" * 50)
|
||||
print("DETAILED SUMMARY STATISTICS")
|
||||
print("=" * 50)
|
||||
print(data['Edep'].describe())
|
||||
```
|
||||
|
||||
Run it and view the figure:
|
||||
|
||||
```bash
|
||||
python3 analyze_advanced.py
|
||||
eog edep_advanced_analysis.png
|
||||
```
|
||||
|
||||
|
||||
|
||||
## Step 7: Comparing Different Particle Types
|
||||
|
||||
Let's generate data for different particles and compare their energy deposition patterns.
|
||||
|
||||
1. Create macro files for different particles. First, protons:
|
||||
|
||||
```bash
|
||||
cd ~/g4-labs/lab3/B1
|
||||
nano proton_beam.mac
|
||||
```
|
||||
|
||||
Add:
|
||||
```
|
||||
/run/numberOfThreads 1
|
||||
/run/initialize
|
||||
/gun/particle proton
|
||||
/gun/energy 150 MeV
|
||||
/run/printProgress 1000
|
||||
/run/beamOn 10000
|
||||
```
|
||||
|
||||
2. Create electron beam macro:
|
||||
|
||||
```bash
|
||||
nano electron_beam.mac
|
||||
```
|
||||
|
||||
Add:
|
||||
```
|
||||
/run/numberOfThreads 1
|
||||
/run/initialize
|
||||
/gun/particle e-
|
||||
/gun/energy 50 MeV
|
||||
/run/printProgress 1000
|
||||
/run/beamOn 10000
|
||||
```
|
||||
|
||||
3. 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.
|
||||
|
||||
```bash
|
||||
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:
|
||||
|
||||
```bash
|
||||
nano compare_particles.py
|
||||
```
|
||||
|
||||
Paste the following content (must start at column 0 — no leading whitespace):
|
||||
|
||||
```python
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# All three CSVs were produced by the modified B1, so they share the
|
||||
# same Geant4 header format. We supply the column name ourselves.
|
||||
read_kwargs = dict(comment='#', header=None, names=['Edep'])
|
||||
|
||||
print("Loading particle data...")
|
||||
gamma_data = pd.read_csv('gamma_6MeV.csv', **read_kwargs)
|
||||
proton_data = pd.read_csv('proton_150MeV.csv', **read_kwargs)
|
||||
electron_data = pd.read_csv('electron_50MeV.csv', **read_kwargs)
|
||||
|
||||
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
|
||||
fig.suptitle('Particle Type Comparison: Energy Deposition',
|
||||
fontsize=16, fontweight='bold')
|
||||
|
||||
ax1 = axes[0, 0]
|
||||
ax1.hist(gamma_data['Edep'], bins=80, alpha=0.7,
|
||||
color='blue', label='6 MeV gamma', edgecolor='black')
|
||||
ax1.set_xlabel('Energy Deposited (MeV)')
|
||||
ax1.set_ylabel('Frequency')
|
||||
ax1.set_title('Gamma Rays (6 MeV)')
|
||||
ax1.legend()
|
||||
ax1.grid(True, alpha=0.3)
|
||||
|
||||
ax2 = axes[0, 1]
|
||||
ax2.hist(proton_data['Edep'], bins=80, alpha=0.7,
|
||||
color='red', label='150 MeV p+', edgecolor='black')
|
||||
ax2.set_xlabel('Energy Deposited (MeV)')
|
||||
ax2.set_ylabel('Frequency')
|
||||
ax2.set_title('Protons (150 MeV)')
|
||||
ax2.legend()
|
||||
ax2.grid(True, alpha=0.3)
|
||||
|
||||
ax3 = axes[1, 0]
|
||||
ax3.hist(electron_data['Edep'], bins=80, alpha=0.7,
|
||||
color='green', label='50 MeV e-', edgecolor='black')
|
||||
ax3.set_xlabel('Energy Deposited (MeV)')
|
||||
ax3.set_ylabel('Frequency')
|
||||
ax3.set_title('Electrons (50 MeV)')
|
||||
ax3.legend()
|
||||
ax3.grid(True, alpha=0.3)
|
||||
|
||||
# Normalized overlay
|
||||
ax4 = axes[1, 1]
|
||||
ax4.hist(gamma_data['Edep'], bins=60, alpha=0.5, color='blue',
|
||||
label=f'gamma (mu={gamma_data["Edep"].mean():.2f})', density=True)
|
||||
ax4.hist(proton_data['Edep'], bins=60, alpha=0.5, color='red',
|
||||
label=f'p+ (mu={proton_data["Edep"].mean():.2f})', density=True)
|
||||
ax4.hist(electron_data['Edep'], bins=60, alpha=0.5, color='green',
|
||||
label=f'e- (mu={electron_data["Edep"].mean():.2f})', density=True)
|
||||
ax4.set_xlabel('Energy Deposited (MeV)')
|
||||
ax4.set_ylabel('Probability Density')
|
||||
ax4.set_title('Normalized Overlay Comparison')
|
||||
ax4.legend()
|
||||
ax4.grid(True, alpha=0.3)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig('particle_comparison.png', dpi=300)
|
||||
print("Comparison figure saved to particle_comparison.png")
|
||||
|
||||
print("\n" + "=" * 60)
|
||||
print("STATISTICAL COMPARISON BY PARTICLE TYPE")
|
||||
print("=" * 60)
|
||||
|
||||
for name, df in [('Gamma (6 MeV)', gamma_data),
|
||||
('Proton (150 MeV)', proton_data),
|
||||
('Electron (50 MeV)', electron_data)]:
|
||||
n_zero = (df['Edep'] == 0).sum()
|
||||
print(f"\n{name}:")
|
||||
print(f" Mean: {df['Edep'].mean():.4f} MeV")
|
||||
print(f" Median: {df['Edep'].median():.4f} MeV")
|
||||
print(f" Std Dev: {df['Edep'].std():.4f} MeV")
|
||||
print(f" Max: {df['Edep'].max():.4f} MeV")
|
||||
print(f" Zero events: {n_zero} ({n_zero/len(df)*100:.1f}%)")
|
||||
```
|
||||
|
||||
Run it:
|
||||
|
||||
```bash
|
||||
python3 compare_particles.py
|
||||
```
|
||||
|
||||
|
||||
|
||||
## Analysis Questions
|
||||
|
||||
Now that you have generated and analyzed the data, answer these questions:
|
||||
|
||||
### Part 1: Understanding the Distribution
|
||||
|
||||
1. Describe the shape of the energy deposition histogram for gamma rays. Is the energy deposition the same for every event?
|
||||
|
||||
2. 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
|
||||
|
||||
3. 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)?
|
||||
|
||||
4. 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
|
||||
|
||||
5. What percentage of events have zero energy deposition? What does this represent physically?
|
||||
|
||||
6. Why would a 6 MeV gamma ray pass through the detector without depositing any energy?
|
||||
|
||||
### Part 3: Comparing Particle Types
|
||||
|
||||
7. Compare the mean energy deposition for:
|
||||
- 6 MeV gammas
|
||||
- 150 MeV protons
|
||||
- 50 MeV electrons
|
||||
|
||||
8. Which particle type deposits the most energy on average? Why?
|
||||
|
||||
9. 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.
|
||||
|
||||
10. 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
|
||||
|
||||
11. In the advanced analysis script, we filtered events with `data[data['Edep'] > 4.0]`. Explain in plain English what this pandas operation does.
|
||||
|
||||
12. What is the advantage of using pandas DataFrames versus manipulating the CSV file directly with basic Python?
|
||||
|
||||
13. 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
|
||||
|
||||
14. Why is CSV a good intermediate format for this analysis workflow (Geant4 → CSV → Python)?
|
||||
|
||||
15. 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:
|
||||
|
||||
- [ ] Modify Geant4 C++ source code to add data output
|
||||
- [ ] Use G4AnalysisManager to create and write CSV files
|
||||
- [ ] Compile modified Geant4 applications
|
||||
- [ ] Generate simulation data with multiple particle types
|
||||
- [ ] Load CSV data into pandas DataFrames
|
||||
- [ ] Calculate statistical measures (mean, std, median, percentiles)
|
||||
- [ ] Create histograms with matplotlib
|
||||
- [ ] Filter and manipulate data using pandas operations
|
||||
- [ ] Create multi-panel figures for comparative analysis
|
||||
- [ ] Interpret energy deposition distributions physically
|
||||
|
||||
## 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.
|
||||
|
||||
## References
|
||||
|
||||
- Geant4 Collaboration. (n.d.). *Geant4 Analysis Tools Documentation*
|
||||
- Hrivnacova, I. (2018). *Using Geant4 Analysis Manager*
|
||||
- Hrivnacova, I. (2023). *Geant4 Analysis Tools Tutorial*
|
||||
- McKinney, W. (2010). *Data Structures for Statistical Computing in Python*
|
||||
- Hunter, J. D. (2007). *Matplotlib: A 2D Graphics Environment*
|
||||
- The Python Graph Gallery. (n.d.). *Histogram Examples*
|
||||
|
||||
## Additional Resources
|
||||
|
||||
- [Pandas Documentation](https://pandas.pydata.org/docs/)
|
||||
- [Matplotlib Gallery](https://matplotlib.org/stable/gallery/index.html)
|
||||
- [Geant4 Analysis Manager Examples](https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Analysis/analysis.html)
|
||||
- [Python Data Science Handbook](https://jakevdp.github.io/PythonDataScienceHandbook/)
|
||||
|
||||
## Optional Extensions
|
||||
|
||||
If you finish early or want to explore further:
|
||||
|
||||
1. **Add More Columns**: Modify the code to save additional data:
|
||||
- Track length
|
||||
- Number of steps
|
||||
- Particle position at end of event
|
||||
|
||||
2. **Different Geometries**: Modify the detector geometry in `DetectorConstruction.cc` and observe how it affects energy deposition
|
||||
|
||||
3. **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
|
||||
|
||||
4. **Interactive Plotting**: Modify the Python scripts to use `plt.show()` for interactive exploration of the data
|
||||
|
||||
5. **Statistical Tests**: Use scipy to perform statistical tests comparing different datasets (t-tests, KS tests, etc.)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 124 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 96 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 399 KiB |
@@ -0,0 +1,138 @@
|
||||
"""Generate annihilation_diagram.png for Lab 4.
|
||||
|
||||
A schematic of positron-electron annihilation in PET:
|
||||
F-18 nucleus -> emits a positron (e+) -> positron travels a short
|
||||
distance in tissue -> annihilates with an electron (e-) -> two 511 keV
|
||||
photons fly off back-to-back.
|
||||
|
||||
Layout: F-18 sits up and to the left, the positron path arcs down to an
|
||||
annihilation point in the middle, and the two 511 keV photons fly off
|
||||
horizontally left and right so they never overlap the positron path.
|
||||
"""
|
||||
|
||||
import os
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from matplotlib.patches import Circle, FancyArrowPatch
|
||||
|
||||
OUT_DIR = os.path.dirname(os.path.abspath(__file__))
|
||||
OUT_PATH = os.path.join(OUT_DIR, "annihilation_diagram.png")
|
||||
|
||||
|
||||
def main() -> None:
|
||||
fig, ax = plt.subplots(figsize=(12, 6.4))
|
||||
ax.set_xlim(0, 12)
|
||||
ax.set_ylim(0, 6.4)
|
||||
ax.set_aspect("equal")
|
||||
ax.axis("off")
|
||||
|
||||
# --- F-18 nucleus, upper left ----------------------------------------
|
||||
nucleus_xy = (1.7, 5.0)
|
||||
ax.add_patch(Circle(nucleus_xy, 0.45, facecolor="#3b6fb6",
|
||||
edgecolor="black", linewidth=1.5, zorder=3))
|
||||
ax.text(*nucleus_xy, "F-18", ha="center", va="center",
|
||||
color="white", fontsize=12, fontweight="bold", zorder=4)
|
||||
ax.text(nucleus_xy[0], nucleus_xy[1] - 0.85,
|
||||
r"$\beta^+$ decay", ha="center", va="center", fontsize=10,
|
||||
style="italic", color="#222")
|
||||
|
||||
# --- Annihilation point, middle-lower --------------------------------
|
||||
annihilation_xy = (6.0, 3.0)
|
||||
|
||||
# --- Wavy positron path from nucleus down to annihilation point ------
|
||||
# Build a parametric curve from F-18 to the annihilation point with a
|
||||
# small sinusoidal wobble that follows the curve direction.
|
||||
t = np.linspace(0.05, 0.97, 240)
|
||||
base_x = nucleus_xy[0] + (annihilation_xy[0] - nucleus_xy[0]) * t
|
||||
base_y = nucleus_xy[1] + (annihilation_xy[1] - nucleus_xy[1]) * t
|
||||
# Tangent direction along the line.
|
||||
dx = annihilation_xy[0] - nucleus_xy[0]
|
||||
dy = annihilation_xy[1] - nucleus_xy[1]
|
||||
seg_len = np.hypot(dx, dy)
|
||||
nx, ny = -dy / seg_len, dx / seg_len # unit normal
|
||||
wobble = 0.18 * np.sin(t * 9.0 * np.pi) * (1.0 - 0.4 * t)
|
||||
path_x = base_x + wobble * nx
|
||||
path_y = base_y + wobble * ny
|
||||
ax.plot(path_x, path_y, color="#d04545", linewidth=1.6, alpha=0.9, zorder=2)
|
||||
ax.annotate(
|
||||
"", xy=annihilation_xy, xytext=(path_x[-3], path_y[-3]),
|
||||
arrowprops=dict(arrowstyle="-|>", color="#d04545", lw=1.6),
|
||||
zorder=2,
|
||||
)
|
||||
|
||||
# Positron label, placed to the side of the path so it doesn't collide
|
||||
# with the photon labels (which sit at y ~= 3.5).
|
||||
midpoint = (path_x[len(path_x) // 2] + 0.45, path_y[len(path_y) // 2] + 0.25)
|
||||
ax.text(*midpoint, r"$e^+$", color="#d04545", fontsize=15,
|
||||
fontweight="bold", ha="center")
|
||||
ax.text(midpoint[0] + 0.6, midpoint[1] - 0.55,
|
||||
"positron range\n(~1 mm in tissue)",
|
||||
ha="left", va="center", fontsize=9, color="#444")
|
||||
|
||||
# --- Annihilation event marker + label -------------------------------
|
||||
ax.plot(*annihilation_xy, marker="*", markersize=46,
|
||||
markerfacecolor="gold", markeredgecolor="black",
|
||||
markeredgewidth=1.4, zorder=3)
|
||||
ax.text(annihilation_xy[0], annihilation_xy[1] - 0.95,
|
||||
"annihilation", ha="center", va="center",
|
||||
fontsize=10, fontweight="bold", color="#222")
|
||||
ax.text(annihilation_xy[0], annihilation_xy[1] - 1.32,
|
||||
r"$e^+ + e^- \to \gamma + \gamma$",
|
||||
ha="center", va="center", fontsize=11, color="#222")
|
||||
|
||||
# --- Inbound electron from below the annihilation point --------------
|
||||
e_minus_xy = (annihilation_xy[0] + 0.95, annihilation_xy[1] - 0.05)
|
||||
ax.text(e_minus_xy[0] + 0.55, e_minus_xy[1] + 0.05,
|
||||
r"$e^-$", color="#1f6f3b", fontsize=15,
|
||||
fontweight="bold", ha="center", va="center")
|
||||
ax.annotate(
|
||||
"", xy=(annihilation_xy[0] + 0.32, annihilation_xy[1]),
|
||||
xytext=(e_minus_xy[0] + 0.30, e_minus_xy[1]),
|
||||
arrowprops=dict(arrowstyle="-|>", color="#1f6f3b", lw=1.6),
|
||||
)
|
||||
|
||||
# --- Two 511 keV photons, horizontal, back-to-back -------------------
|
||||
photon_y = annihilation_xy[1]
|
||||
left_end = (0.6, photon_y)
|
||||
right_end = (11.4, photon_y)
|
||||
|
||||
ax.add_patch(FancyArrowPatch(
|
||||
(annihilation_xy[0] - 0.35, photon_y), left_end,
|
||||
arrowstyle="-|>", mutation_scale=22, color="#b8860b", linewidth=2.2,
|
||||
))
|
||||
ax.add_patch(FancyArrowPatch(
|
||||
(annihilation_xy[0] + 0.35, photon_y), right_end,
|
||||
arrowstyle="-|>", mutation_scale=22, color="#b8860b", linewidth=2.2,
|
||||
))
|
||||
|
||||
# Wavy decoration so the arrows read as electromagnetic waves.
|
||||
for x_lo, x_hi in [(0.8, 5.55), (6.45, 11.2)]:
|
||||
xs = np.linspace(x_lo, x_hi, 240)
|
||||
ys = photon_y + 0.13 * np.sin((xs - x_lo) * 4.5)
|
||||
ax.plot(xs, ys, color="#b8860b", linewidth=1.0, alpha=0.7)
|
||||
|
||||
ax.text(2.8, photon_y + 0.55, r"$\gamma$ 511 keV", color="#8a6308",
|
||||
fontsize=12, fontweight="bold", ha="center")
|
||||
ax.text(9.0, photon_y + 0.55, r"$\gamma$ 511 keV", color="#8a6308",
|
||||
fontsize=12, fontweight="bold", ha="center")
|
||||
|
||||
# --- Footer caption --------------------------------------------------
|
||||
ax.text(6.0, 0.55,
|
||||
"Two 511 keV photons fly back-to-back (180°) to conserve momentum.\n"
|
||||
"A PET ring detects coincident pairs and traces the line connecting "
|
||||
"them — the line of response.",
|
||||
ha="center", va="center", fontsize=10, color="#333")
|
||||
|
||||
ax.set_title(
|
||||
"Positron-Electron Annihilation in PET",
|
||||
fontsize=14, fontweight="bold", pad=12,
|
||||
)
|
||||
|
||||
plt.tight_layout()
|
||||
fig.savefig(OUT_PATH, dpi=200, bbox_inches="tight")
|
||||
print(f"Wrote {OUT_PATH}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 72 KiB |
@@ -0,0 +1,433 @@
|
||||
---
|
||||
title: "Lab 4: Medical Physics Application - A PET Scanner"
|
||||
author: "John Murphy and Steven Carlson"
|
||||
date: "2025-11-11"
|
||||
---
|
||||
|
||||
## 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.
|
||||
|
||||
|
||||
|
||||
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:
|
||||
|
||||
1. Record the energy deposited in each scintillator crystal for every event.
|
||||
2. Plot an energy spectrum and identify the **511 keV photopeak** — the signature of a gamma ray fully absorbed by a single crystal.
|
||||
3. 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
|
||||
|
||||
1. What is positron-electron annihilation, and why does it produce two photons rather than one?
|
||||
2. What is the energy of each annihilation photon, and why is this specific energy important for PET?
|
||||
3. In a detector energy spectrum, what does the "photopeak" represent? What is the "Compton continuum"?
|
||||
4. 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`.
|
||||
|
||||
1. Create a lab directory and copy the B3 source tree into it:
|
||||
|
||||
```bash
|
||||
mkdir -p ~/g4-labs/lab4
|
||||
cd ~/g4-labs/lab4
|
||||
cp -r ~/geant4-v11.3.2/examples/basic/B3 .
|
||||
cd B3/B3a
|
||||
```
|
||||
|
||||
2. Create a build directory and configure with CMake:
|
||||
|
||||
```bash
|
||||
mkdir build
|
||||
cd build
|
||||
cmake ..
|
||||
```
|
||||
|
||||
3. Compile the application:
|
||||
|
||||
```bash
|
||||
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.
|
||||
|
||||
1. Run the application in interactive mode:
|
||||
|
||||
```bash
|
||||
./exampleB3a
|
||||
```
|
||||
|
||||
2. The Qt GUI will open. The geometry is more complex than B1: rings of scintillator crystals surround a central cylinder representing a patient phantom.
|
||||
|
||||
|
||||
|
||||
3. Fire a few events to see particles produced by F-18 decay:
|
||||
|
||||
```bash
|
||||
/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.
|
||||
|
||||
4. Rotate the view to confirm the ring geometry, then exit the GUI:
|
||||
|
||||
```bash
|
||||
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
|
||||
|
||||
1. Open `src/RunAction.cc`:
|
||||
|
||||
```bash
|
||||
gedit src/RunAction.cc &
|
||||
```
|
||||
|
||||
2. Add the analysis-manager header at the top, near the other `#include` lines:
|
||||
|
||||
```cpp
|
||||
#include "G4AnalysisManager.hh"
|
||||
```
|
||||
|
||||
3. Inside `RunAction::BeginOfRunAction(const G4Run*)`, add at the top of the function body:
|
||||
|
||||
```cpp
|
||||
// --- Lab 4 addition: route the framework's score-ntuple output to CSV
|
||||
auto analysisManager = G4AnalysisManager::Instance();
|
||||
analysisManager->SetDefaultFileType("csv");
|
||||
analysisManager->OpenFile("B3_output");
|
||||
```
|
||||
|
||||
4. Inside `RunAction::EndOfRunAction(const G4Run*)`, add at the **very end** of the function:
|
||||
|
||||
```cpp
|
||||
// --- 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
|
||||
|
||||
1. Re-compile your modified application:
|
||||
|
||||
```bash
|
||||
cd build
|
||||
make -j$(nproc)
|
||||
```
|
||||
|
||||
If you encounter errors, double-check brace placement, semicolons, and that your additions are inside the correct function.
|
||||
|
||||
2. Create a macro file `lab4.mac` in the `B3a` directory (one level up from `build`):
|
||||
|
||||
```bash
|
||||
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
|
||||
```
|
||||
|
||||
3. Run the simulation in batch mode:
|
||||
|
||||
```bash
|
||||
cd build
|
||||
./exampleB3a ../lab4.mac
|
||||
```
|
||||
|
||||
This will take a few minutes depending on your hardware.
|
||||
|
||||
4. Verify the output:
|
||||
|
||||
```bash
|
||||
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`:
|
||||
|
||||
```bash
|
||||
nano analyze_spectrum.py
|
||||
```
|
||||
|
||||
Paste the following content (must start at column 0 — no leading whitespace):
|
||||
|
||||
```python
|
||||
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:
|
||||
|
||||
```bash
|
||||
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).
|
||||
|
||||
|
||||
|
||||
## 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:
|
||||
|
||||
```bash
|
||||
nano analyze_coincidence.py
|
||||
```
|
||||
|
||||
Paste the following content (must start at column 0 — no leading whitespace):
|
||||
|
||||
```python
|
||||
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")
|
||||
```
|
||||
|
||||
Run it:
|
||||
|
||||
```bash
|
||||
python3 analyze_coincidence.py
|
||||
eog 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.
|
||||
|
||||
|
||||
|
||||
## Analysis Questions
|
||||
|
||||
### Energy Spectrum
|
||||
|
||||
1. Open `pet_spectrum.png`. At what energy is the main peak centered? How close is this to the theoretical 511 keV?
|
||||
2. Why is the photopeak so much taller than the rest of the spectrum?
|
||||
3. 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
|
||||
|
||||
4. 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*?
|
||||
5. 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
|
||||
|
||||
6. 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.)
|
||||
7. Where in the image is the densest crossing region? What does this tell you about where the simulated source is located?
|
||||
|
||||
### Workflow
|
||||
|
||||
8. 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)?
|
||||
9. 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:
|
||||
|
||||
- [ ] Compile and run the unmodified B3a example
|
||||
- [ ] Identify the existing `crystal/edep` hits collection in the source code
|
||||
- [ ] Modify `RunAction.cc` to route the framework's score ntuples to CSV
|
||||
- [ ] Generate `B3_output_nt_crystal_edep_t0.csv` from a 50,000-event run
|
||||
- [ ] Plot the crystal energy spectrum and identify the 511 keV photopeak
|
||||
- [ ] Map crystal copy numbers to (x, y) positions
|
||||
- [ ] Find true-coincidence events and draw their lines of response
|
||||
|
||||
## 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.
|
||||
- Geant4 Collaboration. (n.d.). *Geant4 Application Developer's Guide — Hits and Sensitive Detectors*.
|
||||
- Cherry, S. R., Sorenson, J. A., & Phelps, M. E. (2012). *Physics in Nuclear Medicine* (4th ed.). Elsevier Health Sciences.
|
||||
|
||||
## Optional Extensions
|
||||
|
||||
1. **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.
|
||||
2. **3-D LORs.** Use the `Z` column from `crystal_position` to draw lines of response in 3-D with `mpl_toolkits.mplot3d`.
|
||||
3. **Two sources.** Modify the primary generator to alternate between two source positions, then check whether your LOR plot shows two overlapping clusters of crossings.
|
||||
4. **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.
|
||||
Reference in New Issue
Block a user