Saving PET Events to ROOT

Table of Contents

Why Save PET Events to ROOT

In a PET simulation, you typically produce many events that contain multiple detector hits, coincidence pairs, and reconstructed lines of response. Looking at text output is not practical for real analysis. The ROOT framework provides efficient storage, fast selection, and powerful plotting for large datasets. Saving PET events from Geant4 to a ROOT file lets you run many different analyses without rerunning the simulation.

In this chapter we focus on how to organize PET data, how to connect it to Geant4’s analysis system, and how to fill ROOT ntuples with the relevant PET quantities. Details of ROOT analysis itself are discussed elsewhere, so here we only prepare the output that ROOT will later read.

Structuring PET Event Data

Before you write any code, decide what information each PET event should contain. For a basic PET scanner example, you usually work at least with coincidences and possibly with the underlying single hits. A typical minimal set of quantities per coincidence is:

Event-level identifiers, such as the Geant4 event ID, to relate coincidences back to Geant4 events.

Coincidence-level information, such as the IDs of the two detector crystals that fired, and the times and energies of the two hits.

Geometrical information, such as the positions of the two detector crystals, or directly the line-of-response (LOR) coordinates.

You can represent each detected hit as a simple C++ structure or class in memory, and then combine two hits into a coincidence structure. For many PET simulations, it is sufficient to store one line per coincidence in the output ROOT ntuple, with columns for all fields that you want to analyze later.

Once you know which fields you want, you must decide where in the Geant4 flow you will create and fill these data structures. Coincidences are properties of an event or of a run, not of an individual tracking step, so the logic for pairing hits is usually located in the event action class or in a dedicated PET manager class that is called from the event action.

Using G4AnalysisManager for PET

Geant4 provides the G4AnalysisManager class to write data in various formats, including ROOT. You already encountered G4AnalysisManager when learning about histograms and ntuples. For the PET example, we will use it to create an ntuple with one row per coincidence.

The common pattern is:

Create and configure G4AnalysisManager in your RunAction constructor, and create your ntuples and columns there.

Open the output file at the beginning of the run, in BeginOfRunAction.

During the simulation, fill histograms and ntuples whenever you have a complete set of values, for example when you identify a coincidence.

At the end of the run, in EndOfRunAction, write the data to disk and close the file.

The important point is that you use the same G4AnalysisManager instance across your user action classes. Typically you retrieve the singleton with:

cpp

auto analysisManager = G4AnalysisManager::Instance();

whenever you need to fill something. You do not create multiple analysis managers.

Important rule: Use a single G4AnalysisManager instance obtained with G4AnalysisManager::Instance(), configure ntuples and histograms in RunAction, open and close the ROOT file in the run actions, and only call FillNtuple*Column and AddNtupleRow after all required columns for that row have been filled.

Defining a PET Ntuple

For PET output, a ROOT ntuple is usually a better choice than only histograms, because you want to apply different cuts and combinations later. You define the ntuple and its columns in the RunAction constructor, after you create the G4AnalysisManager.

A simple example for coincidence data might look like this:

cpp

MyRunAction::MyRunAction()
{
  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->SetVerboseLevel(1);
  analysisManager->SetDefaultFileType("root");
  analysisManager->SetFileName("pet_output");
  // Ntuple ID 0, name "PET", title "PET coincidences"
  analysisManager->CreateNtuple("PET", "PET coincidences");
  // Event and coincidence IDs
  analysisManager->CreateNtupleIColumn("eventID");
  analysisManager->CreateNtupleIColumn("coinID");
  // Detector IDs
  analysisManager->CreateNtupleIColumn("detID1");
  analysisManager->CreateNtupleIColumn("detID2");
  // Energies (in keV, for example)
  analysisManager->CreateNtupleDColumn("energy1");
  analysisManager->CreateNtupleDColumn("energy2");
  // Times (global times, in ns)
  analysisManager->CreateNtupleDColumn("time1");
  analysisManager->CreateNtupleDColumn("time2");
  // Positions of the two hits (world coordinates)
  analysisManager->CreateNtupleDColumn("x1");
  analysisManager->CreateNtupleDColumn("y1");
  analysisManager->CreateNtupleDColumn("z1");
  analysisManager->CreateNtupleDColumn("x2");
  analysisManager->CreateNtupleDColumn("y2");
  analysisManager->CreateNtupleDColumn("z2");
  // Optionally, precomputed LOR direction components or midpoints
  analysisManager->CreateNtupleDColumn("lor_dx");
  analysisManager->CreateNtupleDColumn("lor_dy");
  analysisManager->CreateNtupleDColumn("lor_dz");
  analysisManager->FinishNtuple();
}

This defines a single ntuple with integer and double columns. The detailed set of columns is up to your project. For example, you might add flags for energy window acceptance, ring and block indices instead of raw crystal IDs, or the annihilation point if you keep truth information for validation.

Compared with histograms, an ntuple gives you full flexibility later. You can, for instance, plot an energy spectrum of one detector by selecting detID1 == some_value in ROOT, or create 2D plots of hit position versus energy.

Opening and Closing the ROOT File

The G4AnalysisManager handles writing a ROOT file for you. You control file opening and closing in your run actions, typically once per run. In BeginOfRunAction, you open the file:

cpp

void MyRunAction::BeginOfRunAction(const G4Run*)
{
  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->OpenFile();
}

If you provided a base file name in the constructor, for example "pet_output", this will produce a file such as pet_output.root. You can also choose to include run numbers or other information in the file name, by calling SetFileName before OpenFile if needed.

At the end of the run, you must write and close the file:

cpp

void MyRunAction::EndOfRunAction(const G4Run*)
{
  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->Write();
  analysisManager->CloseFile();
}

If you forget to call Write or CloseFile, the ROOT file may be incomplete or empty. For a PET simulation that can generate many events, it is also a good idea to keep the verbosity level at 0 or 1 to avoid large amounts of terminal output.

Important rule: Always call OpenFile once at the beginning of the run, and both Write and CloseFile at the end of the run. If you do not explicitly write and close, your ROOT file may not contain any data.

Filling Ntuples with PET Coincidences

Once you have a defined ntuple and an open output file, you must decide where to fill rows. For PET, you usually build coincidences from hits recorded in your sensitive detector. The steps look like this:

During an event, sensitive detectors create hits when energy is deposited in crystals. These hits may be stored in a hits collection or directly in a per-event container that you manage.
At the end of the event, when all hits for that event are available, you run coincidence logic. You select valid hits, apply energy windows, and pair hits that are within a coincidence time window.
For each coincidence pair you create, you compute any derived quantities such as energy sums, LOR direction, or midpoints, and then fill one row in the ntuple.

Assume you have a simple coincidence structure:

cpp

struct Coincidence
{
  G4int  eventID;
  G4int  coinID;
  G4int  detID1, detID2;
  G4double energy1, energy2;
  G4double time1, time2;
  G4ThreeVector pos1, pos2;
};

In your EventAction::EndOfEventAction, after you have filled a std::vector<Coincidence> with all coincidences for the event, you can write:

cpp

void MyEventAction::EndOfEventAction(const G4Event* event)
{
  auto analysisManager = G4AnalysisManager::Instance();
  for (const auto& c : fCoincidences) {
    // Optionally compute LOR direction from positions
    G4ThreeVector lorDir = (c.pos2 - c.pos1).unit();
    analysisManager->FillNtupleIColumn(0, c.eventID);
    analysisManager->FillNtupleIColumn(1, c.coinID);
    analysisManager->FillNtupleIColumn(2, c.detID1);
    analysisManager->FillNtupleIColumn(3, c.detID2);
    analysisManager->FillNtupleDColumn(4, c.energy1);
    analysisManager->FillNtupleDColumn(5, c.energy2);
    analysisManager->FillNtupleDColumn(6, c.time1);
    analysisManager->FillNtupleDColumn(7, c.time2);
    analysisManager->FillNtupleDColumn(8,  c.pos1.x());
    analysisManager->FillNtupleDColumn(9,  c.pos1.y());
    analysisManager->FillNtupleDColumn(10, c.pos1.z());
    analysisManager->FillNtupleDColumn(11, c.pos2.x());
    analysisManager->FillNtupleDColumn(12, c.pos2.y());
    analysisManager->FillNtupleDColumn(13, c.pos2.z());
    analysisManager->FillNtupleDColumn(14, lorDir.x());
    analysisManager->FillNtupleDColumn(15, lorDir.y());
    analysisManager->FillNtupleDColumn(16, lorDir.z());
    analysisManager->AddNtupleRow();
  }
  fCoincidences.clear();
}

Column indices must match the creation order in your RunAction. If you change the column order in the ntuple definition, you must update the indices here as well.

In many PET projects, you also want to keep the Geant4 event ID, which you can obtain with:

cpp

G4int eventID = event->GetEventID();

You can store it in your coincidence structures during pairing, so that you can later trace a coincidence back to the original event if needed.

Important rule: For each coincidence, call FillNtuple*Column for every defined column in the same order as the column creation, then finish with AddNtupleRow. Forgetting to call AddNtupleRow means that no data for that coincidence is written to the ntuple.

Choosing What to Store for PET Analysis

You are free to decide how much information to store. More information leads to larger ROOT files but allows more flexible analysis. Less information saves space but may limit what you can study later.

Common choices in a PET example include:

Minimal data per coincidence: detector IDs, energies, times, and the two hit positions. With this, you can create energy spectra, time-of-flight plots, and LORs.

Extended data per coincidence: ring and crystal indices, block indices, acceptance flags from energy or timing windows, and optionally the annihilation point if you have access to generator-level truth. This is very helpful for validation and for debugging the reconstruction.

Full hit lists: instead of only writing coincidences, you can write every single hit in a separate ntuple, or you can add a multiplicity column to your coincidence ntuple and use arrays to store multiple hits. For an absolute beginner course, a simple one-row-per-coincidence ntuple is often the most understandable design.

Table 1 shows an example mapping of PET quantities to ntuple columns:

Concept	Example column name	Type
Geant4 event ID	`eventID`	Integer
Coincidence index	`coinID`	Integer
First detector ID	`detID1`	Integer
Second detector ID	`detID2`	Integer
Energy in first crystal	`energy1`	Double
Energy in second crystal	`energy2`	Double
Time of first hit	`time1`	Double
Time of second hit	`time2`	Double
Position of first hit	`x1, y1, z1`	Double
Position of second hit	`x2, y2, z2`	Double
LOR direction components	`lor_dx, lor_dy, lor_dz`	Double

When you later read the ROOT file, you can easily project any of these columns into histograms and apply cuts such as energy windows or time windows without changing or rerunning the Geant4 simulation.

Working with Multithreading

If your PET simulation uses multithreading, Geant4 also supports parallel output through G4AnalysisManager. In multithreaded mode, each worker thread typically writes its own partial data, and Geant4 merges the outputs when you call Write in the master thread at the end of the run.

At the level of your user code for PET, filling ntuples looks the same in single-threaded and multithreaded runs, as long as you use G4AnalysisManager correctly. You should avoid sharing non-thread-safe containers of coincidences between threads. Instead, each event and thread should handle its own coincidences, and you fill the ntuple directly from EndOfEventAction for that thread.

For simple example applications, you can initially run in single-threaded mode to avoid the additional complexity. Once the logic for creating and filling coincidences is correct, you can enable multithreading and rely on Geant4 to manage most of the details of merging ntuple data.

Summary

Saving PET events to ROOT connects your Geant4 simulation to the analysis work that follows. You decide which PET quantities you want to store, define a ROOT ntuple through G4AnalysisManager, open a ROOT file in the run action, and, for each coincidence, fill one row in the ntuple with detector IDs, energies, times, and positions. With this output, you can use ROOT later to build energy spectra, timing distributions, and image reconstruction studies without modifying or rerunning the simulation.