25.10. Recording Detector Positions

Table of Contents

Why Detector Positions Matter in PET

In a PET simulation you do not only care that a detector fired, you also need to know where it is in space. Each hit must be associated with a detector element, and each detector element has a known position in the scanner ring. From these positions you later build lines of response and, eventually, an image.

In Geant4, the position information exists at two levels that you will typically combine:

The geometrical position of each detector crystal, defined in your geometry (for example in DetectorConstruction).
The interaction position from the simulation step, given by the tracking system when a gamma interacts in the crystal.

For most PET analysis you need both: the crystal index or ID and its nominal center position, and optionally the precise interaction point inside the crystal volume.

In PET reconstruction you must always know which detector element registered the hit and where that element is located in global coordinates. Never record only local positions without a way to convert them to a common global frame.

Geometry, Copy Numbers, and Coordinates

When you build a PET detector ring you typically define one or more volumes repeatedly placed around the patient. Common patterns are a ring of identical modules, each containing an array of crystals. Geant4 lets you distinguish these repeated copies through copy numbers and lets you access their positions in either local or global coordinates.

A PET crystal hit is located by:

The copy number(s) along the geometry hierarchy, which you assign and interpret as detector IDs.
The global position of the step or of the crystal center, which you compute from the geometry transform.

Geant4 uses a right‑handed Cartesian coordinate system. You define detector placements with G4PVPlacement or parameterised volumes. Each placement has a transformation matrix that maps between the local coordinates of a volume and the global coordinates of the world.

The two positions you will most often use are:

The hit position from tracking: G4StepPoint::GetPosition(). This is already in global coordinates.
The geometrical center of a detector element: obtained from its placement transformation, also expressed in global coordinates.

When you talk about the detector position for PET analysis you normally mean the geometrical center of the crystal that fired, in global coordinates. This gives you a stable, reproducible position, independent of the particular interaction point in that crystal.

The position from G4StepPoint::GetPosition() is always in global coordinates. If you want an interaction position in the local coordinates of a crystal you must explicitly transform it yourself.

Storing Detector Positions in Hits

To record detector positions you first need a hits class that can store them. In a PET scanner example you will usually have a custom hit type, for example PetDetectorHit, that extends G4VHit. Along with energy, time, and detector ID, you add a member to store the detector position.

A minimal design might be:

cpp

class PetDetectorHit : public G4VHit {
public:
  PetDetectorHit() = default;
  ~PetDetectorHit() override = default;
  void SetDetectorID(G4int id)      { fDetID = id; }
  void SetPosition(const G4ThreeVector& p) { fPos = p; }
  void SetTime(G4double t)          { fTime = t; }
  void SetEnergy(G4double e)        { fEnergy = e; }
  G4int          GetDetectorID() const { return fDetID; }
  const G4ThreeVector& GetPosition() const { return fPos; }
  G4double       GetTime() const     { return fTime; }
  G4double       GetEnergy() const   { return fEnergy; }
private:
  G4int          fDetID = -1;
  G4ThreeVector  fPos;
  G4double       fTime   = 0.;
  G4double       fEnergy = 0.;
};

Within your sensitive detector, you create a new hit for each relevant energy deposition. You then fill the detector ID and position information in ProcessHits.

You have two natural choices for what to store in fPos:

The interaction point inside the crystal, from the step position.
The crystal center from geometry, independent of the step location.

For PET event building and lines of response, the crystal center is usually more convenient. You can still keep the interaction point separately if you need detailed studies of positron range or depth of interaction.

Always store positions using Geant4 units, for example mm. Do not convert to unit‑less values or external units inside your hits classes. Convert only when you write out or analyze the data.

Getting Crystal Positions from the Geometry

To obtain the nominal position of a detector crystal you need to:

Identify which volume was hit.
Get the touchable history for that step.
Read the copy numbers for your ring and crystal indices.
Compute the global position of the crystal center from the touchable transform.

The touchable encapsulates the full placement information from the world down to the current volume. You access it inside ProcessHits using the pre step point:

cpp

G4StepPoint* prePoint = step->GetPreStepPoint();
const G4TouchableHandle& touchable = prePoint->GetTouchableHandle();

You can now obtain copy numbers along the hierarchy. For example, if your geometry has a ring of modules and each module has an array of crystals, you could assign copy numbers like this during construction:

Level 0, the world volume.
Level 1, the ring or module placement, copy number = module ID.
Level 2, the crystal placement inside the module, copy number = crystal index.

In ProcessHits you can recover them:

cpp

G4int moduleID  = touchable->GetCopyNumber(1);
G4int crystalID = touchable->GetCopyNumber(0);

The indices depend on how many levels you placed below the world, and you must keep your own convention consistent between geometry construction and hit processing.

To get the global position of the crystal center, use the translation part of the transform that places the current volume:

cpp

G4ThreeVector crystalCenter =
  touchable->GetHistory()
           ->GetTopTransform()
           .Inverse()
           .TransformPoint(G4ThreeVector(0., 0., 0.));

Here G4ThreeVector(0,0,0) is the local origin of the crystal volume. Applying the inverse of the top transform converts this local point to global coordinates. The result is the nominal center of the crystal in the world frame.

You then store this position into your hit:

cpp

auto hit = new PetDetectorHit;
hit->SetDetectorID(crystalID);      // or some combined ID that includes module
hit->SetPosition(crystalCenter);

If you also want the precise interaction position you can take it directly from the step:

cpp

G4ThreeVector interactionPos = prePoint->GetPosition(); // already global

You may store both positions in your hit, for instance as two separate members, if you intend to study effects such as depth of interaction.

The GetCopyNumber(level) indices count from the bottom of the geometry hierarchy, not from the world. Always check your volume tree if the IDs look wrong, and keep a clear mapping between copy numbers and your PET detector indices.

Recording Positions in Analysis Output

Once detector positions are stored in hits, you still need to write them to an output format for later analysis, for example a ROOT file. In a PET application you typically do this in EventAction or RunAction, after the tracking is done, by looping over the hit collections and filling an analysis object such as an ntuple.

First, declare suitable ntuple columns in your analysis setup, for example in RunAction:

cpp

auto analysisManager = G4AnalysisManager::Instance();
analysisManager->CreateNtuple("PET", "PET events");
analysisManager->CreateNtupleIColumn("eventID");
analysisManager->CreateNtupleIColumn("detID");
analysisManager->CreateNtupleDColumn("Edep");
analysisManager->CreateNtupleDColumn("time");
analysisManager->CreateNtupleDColumn("x");
analysisManager->CreateNtupleDColumn("y");
analysisManager->CreateNtupleDColumn("z");
analysisManager->FinishNtuple();

Then, at the end of each event, you retrieve the hit collection and fill one ntuple row per hit:

cpp

void EventAction::EndOfEventAction(const G4Event* event)
{
  auto hce = event->GetHCofThisEvent();
  if (!hce) return;
  auto hcID = fHitsCollectionID; // stored earlier from SD registration
  auto hitsCollection =
    static_cast<PetDetectorHitsCollection*>(hce->GetHC(hcID));
  if (!hitsCollection) return;
  auto analysisManager = G4AnalysisManager::Instance();
  G4int eventID = event->GetEventID();
  for (const auto hit : *hitsCollection) {
    G4ThreeVector pos = hit->GetPosition();
    analysisManager->FillNtupleIColumn(0, eventID);
    analysisManager->FillNtupleIColumn(1, hit->GetDetectorID());
    analysisManager->FillNtupleDColumn(2, hit->GetEnergy());
    analysisManager->FillNtupleDColumn(3, hit->GetTime());
    analysisManager->FillNtupleDColumn(4, pos.x() / mm);
    analysisManager->FillNtupleDColumn(5, pos.y() / mm);
    analysisManager->FillNtupleDColumn(6, pos.z() / mm);
    analysisManager->AddNtupleRow();
  }
}

Here the coordinates are converted to millimeters as plain numbers before writing. You can apply any convention you like, but you must be consistent with the units you use in downstream analysis.

For PET, this output provides everything you need to form lines of response: the event identifier, the detector ID, and the global position of the detector element that registered the hit.

Use a clear and consistent unit convention in your output files. Document which units you use for position, time, and energy, and keep them fixed across runs so that downstream PET analysis scripts do not silently misinterpret your data.