25.13 Analyzing PET Data with ROOT

Table of Contents

Preparing to Analyze PET Data

In a PET simulation you typically end with a ROOT file that contains event-level information from Geant4, such as detector IDs, times, energies, and line of response parameters. This chapter focuses on how to work with that ROOT output, not on how to generate it in Geant4. The goal is to go from a collection of coincidence events to physics and imaging quantities that you can inspect, plot, and understand.

Although PET analysis can become very sophisticated, you can already learn a lot from simple ROOT macros that read your trees, apply basic cuts, and fill histograms.

Exploring the ROOT File

Start by opening your PET ROOT file with the ROOT browser. If your Geant4 analysis used G4AnalysisManager to create an ntuple called for example "PET", you will usually find a TTree like "PET" or "PET_ntuple" inside the file.

In a ROOT session, open the file and list its contents:

cpp

TFile *f = TFile::Open("pet_output.root");
f->ls();

Look for a TTree that contains one entry per coincidence event or per detected hit, depending on how you set up your Geant4 output. For PET analysis it is convenient to have a tree where each entry corresponds to one coincidence event, with both hits stored in the same row. A simple structure could be:

Branch name	Type	Meaning
`E1`	Double	Energy in detector 1
`E2`	Double	Energy in detector 2
`t1`	Double	Time of detector 1 hit
`t2`	Double	Time of detector 2 hit
`detID1`	Int	ID of detector 1
`detID2`	Int	ID of detector 2
`x1`, `y1`, `z1`	Double	Position of detector 1 hit
`x2`, `y2`, `z2`	Double	Position of detector 2 hit
`LORx`, `LORy`, `LORz`	Double	Derived line of response parameters

The exact branch names and structure depend on your Geant4 implementation, but the analysis ideas are the same. Use the TBrowser in ROOT to inspect branches interactively:

cpp

TBrowser *b = new TBrowser();

Then double click on your TTree and explore the branches to confirm what is stored.

Basic Quality and Sanity Checks

Before doing detailed analysis, check that the data are sensible. Plot simple 1D and 2D distributions directly with ROOT.

For example, inspect the energy deposited in one crystal:

cpp

TTree *t = (TTree*)f->Get("PET");
t->Draw("E1");

You should see a distribution reflecting the gamma interaction spectrum in a single crystal. If you are simulating 511 keV photons, the spectrum will show a Compton continuum and maybe a visible photopeak if your geometry and materials allow.

Check correlations between the energies in the two detectors:

cpp

t->Draw("E2:E1","","colz");

You expect a band roughly centered around the 511 keV line if you performed coincidence selection already, with events where both energies are around the photopeak and many with one or both energies lower due to partial energy deposition.

Perform similar sanity checks for hit positions and times:

cpp

t->Draw("t2-t1");
t->Draw("z1");
t->Draw("z2");

These quick checks often reveal obvious problems like uninitialized branches, units misunderstandings, or missing cuts in your simulation output.

Always confirm units and coordinate conventions before interpreting PET results. Check that energies are in keV or MeV as expected, times are in ns or ps consistently, and positions are in mm or cm in agreement with your Geant4 geometry.

Energy Windowing in ROOT

PET relies on selecting gamma interactions that are likely to correspond to true 511 keV annihilation photons. This is done with an energy window around the photopeak. Even if you applied an energy window in Geant4, it is useful to be able to change or refine that window in ROOT.

A typical PET energy window might be something like 350 keV to 650 keV around the 511 keV line. In ROOT you can apply this window dynamically using TTree::Draw.

For a symmetric window on both detectors you could write:

cpp

double Emin = 350.0; // keV
double Emax = 650.0; // keV
TString cut = Form("E1>%f && E1<%f && E2>%f && E2<%f", Emin, Emax, Emin, Emax);
t->Draw("E1", cut);

You can then overlay distributions for different window choices to study the effect of tightening or loosening the energy selection. For example, compare two windows:

cpp

t->Draw("E1", "E1>400 && E1<600", "");
t->Draw("E1", "E1>450 && E1<550", "same");

The narrower window will reduce scatter events but also reduce statistics. You can quantify the fraction of events that pass the window by using TTree::GetEntries:

cpp

Long64_t total = t->GetEntries();
Long64_t passed = t->GetEntries(cut);
double efficiency = double(passed) / double(total);

This gives a simple measure of the detection efficiency versus energy window width, which is often interesting when you study trade-offs in PET system design.

Coincidence and Time-of-Flight Analysis

If you stored event-level coincidence information in your ROOT TTree, you can analyze the timing difference between the two detectors in a coincidence pair. This is important both for validating the coincidence selection and for exploring time-of-flight capabilities.

Compute the time difference:

cpp

t->Draw("t2 - t1", cut);

where cut is the energy window condition you defined earlier. The resulting histogram should show a peak centered around zero if your time stamps are correctly referenced. The width of this peak is related to the time resolution of your detector system.

You can fit this peak with a Gaussian:

cpp

TH1F *hDT = new TH1F("hDT","Time difference; t2-t1 [ns]; Counts",200,-5,5);
t->Draw("t2-t1 >> hDT", cut);
hDT->Fit("gaus");

The Gaussian sigma gives a simple estimate of the coincidence timing resolution. If you are simulating time-of-flight PET, you can then relate time difference to position along the line of response. For annihilation photons traveling at speed of light $c$, the position offset $\Delta x$ from the middle of the detector pair satisfies:

$$
\Delta x = \frac{c \Delta t}{2}
$$

where $\Delta t = t_2 - t_1$.

The relation between time difference and position along the line of response is
$$
\Delta x = \frac{c \Delta t}{2}
$$
Make sure your time units and $c$ units are consistent, for example $\Delta t$ in ns and $c \approx 0.3\ \text{mm}/\text{ps} = 300\ \text{mm}/\text{ns}$.

In ROOT you can create a histogram of the inferred position along the LOR from the time difference:

cpp

const double c_mm_per_ns = 299.792458; // mm/ns
t->Draw(Form("0.5*%f*(t2-t1)", c_mm_per_ns), cut);

This type of plot is mainly useful to understand the intrinsic timing spread and how it would translate into spatial uncertainty along the LOR.

Simple Image Reconstruction Plots

A full PET image reconstruction algorithm is beyond the scope of a beginner course. However, you can still create intuitive 2D maps that show where coincidence events are likely to originate in space.

For a simple ring PET geometry with detectors around the origin in the transverse plane, and a small object placed near the center, you can approximate the event origin as the midpoint of the two detector hit positions. This is a rough backprojection method, but it helps visualize the emission distribution.

For each coincidence event compute the midpoint coordinates:

$$
x_{\text{mid}} = \frac{x_1 + x_2}{2},\quad
y_{\text{mid}} = \frac{y_1 + y_2}{2}
$$

You can fill these directly in ROOT:

cpp

TString cut = "E1>350 && E1<650 && E2>350 && E2<650";
t->Draw("(y1+y2)/2.0:(x1+x2)/2.0", cut, "colz");

The resulting 2D histogram shows a simple backprojected activity map. For a point source at the center of the ring you should see a localized region of higher counts around the center. For an extended phantom the pattern will approximate its shape.

You can also restrict the analysis to a specific axial slice if you stored z1 and z2. For example, select only events where both detectors are in a central slice:

cpp

cut += " && abs(z1) < 5 && abs(z2) < 5";
t->Draw("(y1+y2)/2.0:(x1+x2)/2.0", cut, "colz");

This demonstrates the connection between detector hits and reconstructed image space, even if it is not a proper tomographic reconstruction.

Evaluating PET Performance Metrics

With your ROOT data you can estimate basic performance quantities of your simulated PET system. These include simple count-based measures like true, scatter, and random coincident event fractions, and geometric relationships such as sinograms.

If you tagged events in Geant4 with information about whether they are true, scattered, or random coincidences, you can store a classification branch, for example:

Branch	Type	Meaning
`type`	Int	0 true, 1 scatter, 2 random, etc.

Then in ROOT you can count how many events of each type pass your energy and timing cuts:

cpp

TString commonCut = "E1>350 && E1<650 && E2>350 && E2<650";
Long64_t nTrue = t->GetEntries(commonCut + " && type==0");
Long64_t nScatter = t->GetEntries(commonCut + " && type==1");
Long64_t nRandom = t->GetEntries(commonCut + " && type==2");
double total = double(nTrue + nScatter + nRandom);
double fTrue = nTrue / total;
double fScatter = nScatter / total;
double fRandom = nRandom / total;

You can then print and compare these fractions for different choices of energy windows, timing windows, or geometry configurations.

Another important PET concept is the sinogram, which represents the number of events as a function of LOR projection angle and radial offset. If your detectors are arranged in a simple ring and you stored detector IDs that map to angular positions, you can try a basic sinogram creation.

Assume that each detector ID corresponds to angle $\phi_i$ around the ring, and you have a mapping array:

cpp

const int nDet = 64;
double angle[nDet];
for (int i = 0; i < nDet; ++i) {
  angle[i] = 2.0*TMath::Pi()*i/nDet; // simple ring
}

In a ROOT macro, you can loop over entries and fill a 2D histogram for sinogram bins. A very simple geometry might use the angle of the LOR and a crude radial coordinate, although full sinogram definition can be more detailed.

Even without full sinogram reconstruction, counting coincidences as a function of detector pair angle is informative:

cpp

TH1F *hAngle = new TH1F("hAngle","LOR angle; angle [rad]; Counts", 64, 0, 2*TMath::Pi());
int det1, det2;
t->SetBranchAddress("detID1", &det1);
t->SetBranchAddress("detID2", &det2);
Long64_t nEntries = t->GetEntries();
for (Long64_t i = 0; i < nEntries; ++i) {
  t->GetEntry(i);
  double phi1 = angle[det1];
  double phi2 = angle[det2];
  double lorAngle = 0.5*(phi1 + phi2);
  if (lorAngle < 0) lorAngle += 2.0*TMath::Pi();
  if (lorAngle >= 2.0*TMath::Pi()) lorAngle -= 2.0*TMath::Pi();
  hAngle->Fill(lorAngle);
}
hAngle->Draw();

For a symmetric central point source you expect a flat distribution in angle. Deviations from flatness can indicate nonuniform detection efficiency, geometric effects, or simulation issues.

Writing Simple ROOT Macros for PET

Rather than issuing many TTree::Draw commands interactively, it is convenient to encapsulate your analysis steps in simple ROOT macros. This also improves reproducibility and makes it easier to share your analysis with others.

A minimal PET analysis macro might follow this pattern:

cpp

void analyze_pet(const char* filename="pet_output.root") {
  TFile *f = TFile::Open(filename);
  if (!f || f->IsZombie()) {
    printf("Cannot open file %s\n", filename);
    return;
  }
  TTree *t = (TTree*)f->Get("PET");
  if (!t) {
    printf("Cannot find tree PET\n");
    return;
  }
  // Define cuts
  TString cutEnergy = "E1>350 && E1<650 && E2>350 && E2<650";
  // Energy spectra
  TCanvas *c1 = new TCanvas("c1","Energy spectra",800,600);
  t->Draw("E1>>hE1(200,0,700)", cutEnergy);
  t->Draw("E2>>hE2(200,0,700)", cutEnergy, "same");
  // Time difference
  TCanvas *c2 = new TCanvas("c2","Time difference",800,600);
  t->Draw("t2-t1>>hDT(200,-5,5)", cutEnergy);
  hDT->Fit("gaus");
  // Simple backprojection
  TCanvas *c3 = new TCanvas("c3","Backprojection",800,600);
  t->Draw("(y1+y2)/2.0:(x1+x2)/2.0>>hXY(200,-100,100,200,-100,100)",
          cutEnergy, "colz");
  // Save canvases
  c1->SaveAs("pet_energy.png");
  c2->SaveAs("pet_timing.png");
  c3->SaveAs("pet_backprojection.png");
}

Running:

cpp

root -l -q 'analyze_pet.C("pet_output.root")'

creates and saves the key plots in one step. You can then iterate on this macro to add more performance metrics, such as efficiency curves, scatter fractions, or angle distributions.

Keep your PET analysis macros under version control and document the branches and cuts they use. This ensures that you can reproduce figures, track changes, and avoid silent mismatches between simulation versions and analysis code.

By combining these ROOT techniques with the PET-specific information you already recorded in your Geant4 simulation, you gain a flexible environment for exploring detector performance, event selection strategies, and basic imaging properties.