25.11 Creating Lines of Response

Table of Contents

Concept of a Line of Response

In a PET scanner, each valid coincidence event corresponds to the (approximate) detection of two annihilation photons that were emitted back to back. Each photon is detected in a different detector crystal. If you know the positions of these two detector elements, you can construct the straight line that connects them. This line is called a line of response, often abbreviated as LOR.

In your Geant4 simulation, creating LORs means transforming the information you already have for coincidence events into a simple geometric description: for each coincidence, you store the coordinates of the two detector hits and treat the line between them as one measured LOR. Later, reconstruction code can use these LORs to estimate the spatial distribution of activity in the field of view.

For the purpose of this beginner example, an LOR will usually be represented in the simplest possible way, by the two 3D points where the photons were detected. More compact parameterizations are possible, but the two end points are straightforward to compute and easy to write to an output file.

An LOR in this context is defined by two distinct detector hit positions that belong to the same coincidence event and are not located in the same detector element.

Choosing an LOR Representation

The simplest representation of a line of response uses two global position vectors, one for each detector hit. In Geant4, these positions are typically stored as $G4ThreeVector$ objects in your hit class or in your coincidence event class.

A basic representation in C++ might conceptually look like this:

cpp

class PetLOR {
public:
  G4ThreeVector p1;  // first detector hit position, global coordinates
  G4ThreeVector p2;  // second detector hit position, global coordinates
  G4double      t1;  // detection time of first hit
  G4double      t2;  // detection time of second hit
  G4double      e1;  // deposited energy in first detector
  G4double      e2;  // deposited energy in second detector
};

You do not have to create a dedicated C++ class if you prefer to write values directly into an ntuple, but it is often clearer to think in terms of an object like this. The essential part of the LOR is the pair of positions $(\mathbf{p}_1, \mathbf{p}_2)$. Time and energy can be stored with the LOR for later use, for example for time of flight analysis or for applying additional energy cuts offline.

It is useful to emphasize that you do not need to compute a unit direction vector or parametric form of the line unless you specifically need it. The reconstruction software can always compute such derived quantities from the two endpoints:
$$
\mathbf{d} = \mathbf{p}_2 - \mathbf{p}_1, \quad
\hat{\mathbf{u}} = \frac{\mathbf{d}}{|\mathbf{d}|}.
$$

The minimal LOR data that you must store is:

Position of hit 1, $\mathbf{p}_1 = (x_1, y_1, z_1)$.
Position of hit 2, $\mathbf{p}_2 = (x_2, y_2, z_2)$.
Both positions should be in a consistent coordinate system, usually the global coordinates of the simulation.

Extracting Detector Positions from Hits

Before you can form an LOR, you must extract the positions of the two detector hits that form a coincidence. Earlier in the PET example, you recorded detector hits and then selected coincidence events. At the stage where you know which two hits belong to a coincidence, you already have access to each hit’s position, time, energy, and detector ID.

A typical PET hit class will have at least the following members:

cpp

class PetHit : public G4VHit {
public:
  G4ThreeVector fPosition;
  G4double      fTime;
  G4double      fEnergy;
  G4int         fDetectorID;
  // ...
};

When coincidence finding logic has identified two hits, say hitA and hitB, you obtain the positions like this:

cpp

G4ThreeVector p1 = hitA->GetPosition();
G4ThreeVector p2 = hitB->GetPosition();

These positions should already be in global coordinates if they were recorded that way in the hit. If you stored local coordinates instead, you must convert them to global coordinates by applying the relevant transform for the detector volume, but for a beginner PET example it is usually simpler to store global positions from the beginning.

It is also common to read out associated information at the same time:

cpp

G4double t1 = hitA->GetTime();
G4double t2 = hitB->GetTime();
G4double e1 = hitA->GetEnergy();
G4double e2 = hitB->GetEnergy();
G4int detID1 = hitA->GetDetectorID();
G4int detID2 = hitB->GetDetectorID();

You can optionally apply final checks here, for example to ensure that the two hits are not in the same detector crystal, and that they satisfy any time and energy conditions you require.

Always verify that:

The two hits belong to the same PET event,
They passed your energy window and timing criteria,
Their detector IDs are different.
Only then should you construct an LOR.

Building the LOR from Coincidence Hits

Once you have two valid hits, the next step is to build the LOR and prepare it for output. Conceptually, you take the two global positions and store them together as one line of response.

If you use a dedicated LOR object, the creation is straightforward:

cpp

PetLOR lor;
lor.p1 = p1;
lor.p2 = p2;
lor.t1 = t1;
lor.t2 = t2;
lor.e1 = e1;
lor.e2 = e2;

If your focus is on simple PET without time of flight, you can already use this. If you plan to explore time of flight PET, you may also want to compute time and position along the line, but that can be done later in the analysis code.

A useful derived quantity is the midpoint of the LOR, which is often used for simple backprojection methods:

$$
\mathbf{m} = \frac{\mathbf{p}_1 + \mathbf{p}_2}{2}.
$$

In C++ this is:

cpp

G4ThreeVector mid = 0.5 * (p1 + p2);

Storing this midpoint along with the endpoints can simplify very basic reconstruction demonstrations, because you can, for example, fill a 2D histogram of $(x_{\text{mid}}, y_{\text{mid}})$ to visualize where events concentrate.

Another useful derived quantity, especially for a symmetric ring, is the LOR length:

$$
L = |\mathbf{p}_2 - \mathbf{p}_1|.
$$

You can use this to perform sanity checks, for instance by comparing $L$ to the expected diameter of your detector ring.

The LOR direction and midpoint follow directly from the endpoints:
\[
\mathbf{d} = \mathbf{p}_2 - \mathbf{p}_1,\quad
L = |\mathbf{d}|,\quad
\hat{\mathbf{u}} = \frac{\mathbf{d}}{L},\quad
\mathbf{m} = \frac{\mathbf{p}_1 + \mathbf{p}_2}{2}.
\]
You do not need to store all of these if storage space is limited. The endpoints are sufficient.

Storing LORs with G4AnalysisManager

To use the LORs outside Geant4, for example in ROOT-based reconstruction, you must write their information to an output file. In this course, the Geant4 analysis system with G4AnalysisManager is used to record such data. You have already used it for energies and positions. Now you extend your ntuple or histograms to include LOR information.

A convenient way is to create an ntuple dedicated to LORs, with one row per coincidence event. A minimal structure might look like the following table.

Column name	Content
eventID	Event ID in the Geant4 simulation
x1, y1, z1	Coordinates of hit 1
x2, y2, z2	Coordinates of hit 2
t1, t2	Times of hit 1 and hit 2
e1, e2	Deposited energies in each detector

In your analysis initialization you define this ntuple:

cpp

auto analysisManager = G4AnalysisManager::Instance();
analysisManager->CreateNtuple("LOR", "PET lines of response");
analysisManager->CreateNtupleIColumn("eventID");
analysisManager->CreateNtupleDColumn("x1");
analysisManager->CreateNtupleDColumn("y1");
analysisManager->CreateNtupleDColumn("z1");
analysisManager->CreateNtupleDColumn("x2");
analysisManager->CreateNtupleDColumn("y2");
analysisManager->CreateNtupleDColumn("z2");
analysisManager->CreateNtupleDColumn("t1");
analysisManager->CreateNtupleDColumn("t2");
analysisManager->CreateNtupleDColumn("e1");
analysisManager->CreateNtupleDColumn("e2");
analysisManager->FinishNtuple();

Then, when your PET coincidence code has just created an LOR, you fill these columns. This usually takes place in an event or run action, after you have performed coincidence sorting for the event:

cpp

G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
analysisManager->FillNtupleIColumn(0, eventID);
analysisManager->FillNtupleDColumn(1, lor.p1.x());
analysisManager->FillNtupleDColumn(2, lor.p1.y());
analysisManager->FillNtupleDColumn(3, lor.p1.z());
analysisManager->FillNtupleDColumn(4, lor.p2.x());
analysisManager->FillNtupleDColumn(5, lor.p2.y());
analysisManager->FillNtupleDColumn(6, lor.p2.z());
analysisManager->FillNtupleDColumn(7, lor.t1);
analysisManager->FillNtupleDColumn(8, lor.t2);
analysisManager->FillNtupleDColumn(9, lor.e1);
analysisManager->FillNtupleDColumn(10, lor.e2);
analysisManager->AddNtupleRow();

The exact column indices depend on the order in which you defined them, so make sure they are consistent. Once the run finishes, the analysis manager will write these LOR entries to a file format of your choice, for example a ROOT file. Later chapters will show how to open this file with ROOT and perform LOR based reconstruction steps.

When using G4AnalysisManager to record LORs:

Create a dedicated ntuple for LORs.
Call FillNtuple*Column only after you have confirmed a valid coincidence.
Call AddNtupleRow() once for each LOR you create.

Sanity Checks and Simple Diagnostics

Before relying on your LOR data for reconstruction, it is important to confirm that your simulation is writing reasonable information. A few simple diagnostics can be applied directly to the stored LORs.

You can check that the LOR endpoints lie on the detector ring by plotting their radius in ROOT. If your detector ring is centered at the origin and lies in the $x$-$y$ plane, the transverse radius of a point is
$$
r = \sqrt{x^2 + y^2}.
$$
If the crystals are arranged at radius $R$, most endpoints should cluster near $R$.

You can also inspect the distribution of LOR midpoints. For a point-like positron source at the center, the midpoint distribution should cluster near the origin. If you move the source in the simulation, you should see the cluster follow the source position.

Finally, you can check that the LOR directions look isotropic as seen from the source. For example, you can plot the polar angle and azimuthal angle of the vector $\mathbf{d} = \mathbf{p}_2 - \mathbf{p}_1$.

Performing these simple plots on the LOR file is an effective way to verify that your coincidence selection and LOR construction are correct before you proceed to more complex PET image reconstruction.