KAHIBARO
Discord Login Register

16.4. Energy Spectra

Creating energy distributions

In many detector simulations you are not only interested in the total deposited energy per event, but also in how those energies are distributed over many events. This distribution is what we usually call an energy spectrum. In Geant4, building an energy spectrum is conceptually simple: for each event you compute one or more energy values of interest, then you fill those values into histograms. After many events, the shape of the histogram represents the energy distribution.

Energy spectra are typically created in user action classes that already collect energy information, such as EventAction or SteppingAction, and are written out through the Geant4 analysis system.

From deposited energy to spectral bins

To create an energy spectrum, you first need a definition of the quantity you want to plot. Common examples are the total energy deposited in a detector per event, the energy deposited in a particular detector element, or the energy of particles escaping a region.

For a simple detector, you usually accumulate the total deposited energy in a member variable of EventAction. Each step through a sensitive volume contributes an incremental energy deposit, which you obtain from the step object. The typical pattern is to initialize the per event energy to zero at the beginning of the event, add the deposited energy from each relevant step, then use the final sum at the end of the event when filling the spectrum.

In SteppingAction, you can access the deposited energy for a step through the step object using its total energy deposit. The following C++ call returns the deposited energy in Geant4 internal units, most often MeV:

cpp
G4double edep = step->GetTotalEnergyDeposit();

This energy is usually added to an accumulator that lives in EventAction or in a thread local container. At the end of the event, the accumulated value is a single number that is convenient to use as an entry in a histogram.

Important rule: The value returned by GetTotalEnergyDeposit() is in Geant4 internal units. Always convert to the desired output unit, for example edep / MeV, before filling histograms or ntuples, so that the spectrum axis has a clear physical meaning.

Once you have a scalar energy per event, you are ready to build an energy distribution by filling a histogram.

Defining histograms for spectra

Energy spectra in Geant4 are usually represented as one dimensional histograms. These are created through the analysis manager, which hides the details of the underlying output library.

You typically define your histograms during the initialization phase of the run, usually in RunAction. First you obtain the singleton analysis manager instance, then you create a one dimensional histogram with a chosen name, title, number of bins, and energy range. For example, to record an energy spectrum in MeV between 0 and 10 MeV with 1000 bins, you could write:

cpp
auto analysisManager = G4AnalysisManager::Instance();
G4int id = analysisManager->CreateH1(
  "Edep",                           // histogram name
  "Energy deposition in detector",  // histogram title
  1000,                             // number of bins
  0.,                               // lower edge in MeV
  10.                               // upper edge in MeV
);

The histogram identifier returned by CreateH1 is needed later when filling the spectrum. The range and the number of bins determine the resolution and statistical precision of the spectrum. If the range is too small, you will lose higher energy events. If the number of bins is too small, important features such as peaks or edges may be smeared out. On the other hand, too many bins can give very noisy spectra for limited statistics.

A practical way to choose the bin width is to relate it to the detector resolution. For example, if your detector energy resolution is about 100 keV at a given energy, there is little benefit in choosing a bin width much smaller than 10 keV, since the detector response itself is broader.

Important guideline: Choose histogram binning so that the bin width is smaller than, but not much smaller than, the detector resolution. This avoids both oversampling of noisy data and the loss of relevant spectral structure.

Filling histograms with per event energy

Once your histograms are defined, you need to fill them for each event. The usual place to do this is at the end of an event, when the total energy deposition for that event is known.

In EventAction, you might keep a member variable such as fEdep that accumulates the total energy deposit during the event. In the method that is called at the end of the event, you convert this energy to your chosen unit and fill the histogram through the analysis manager:

cpp
void MyEventAction::EndOfEventAction(const G4Event*) {
  auto analysisManager = G4AnalysisManager::Instance();
  // Convert to MeV for the spectrum
  G4double edepMeV = fEdep / MeV;
  analysisManager->FillH1(0, edepMeV);  // 0 is the histogram ID
}

You can create and fill more than one histogram. For instance, you might want a spectrum for each detector element, or separate spectra for different regions of your geometry.

Sometimes you may want finer control, for example filling a spectrum only if a certain condition is met, such as a coincidence requirement or a particular particle type. In such cases, you apply the cuts in your EventAction or other user action and fill the histogram only when the event passes the selection.

Table of typical energy spectra you might create:

Spectrum typeQuantity per eventTypical source class
Total detector energy spectrumSum of all energy deposits in a volumeEventAction
Per crystal energy spectrumEnergy in each detector elementEventAction with IDs
Escaping particle energy spectrumKinetic energy of exiting particlesSteppingAction or TrackingAction
Region specific dose-related spectrumEnergy per unit mass in a specific regionCustom scoring or SteppingAction

For each of these, the same basic principle applies. You identify a scalar value that characterizes the event or track, and then you fill that value into the appropriate histogram.

Multiple spectra and detector elements

In setups with many detector elements, such as arrays of crystals or segmented calorimeters, you often want a separate spectrum for each element. There are two common approaches.

The first approach is to create one histogram per detector element. During initialization you loop over your detector indices and call CreateH1 with a unique name for each crystal. The histogram identifiers are stored, for example in a vector, and during EventAction you fill the appropriate histogram based on the detector ID. This is simple to understand, but the number of histograms can become large.

The second approach is to use ntuples. For each event, you record both the detector ID and the energy, and you fill an ntuple row with these values. After the run, you can use external tools such as ROOT to build spectra offline. For beginners, one histogram per element is usually more intuitive when the number of elements is modest.

In either case, the logic to connect an energy deposition to a particular detector element usually involves obtaining a copy number or a detector ID from the touched logical volume or from a hit object. You should already have this information available if you have implemented sensitive detectors earlier in the course.

Key point: Always keep a clear mapping between the physical detector element and its corresponding spectrum. Use consistent detector IDs and histogram indices to avoid mixing signals from different parts of your detector.

Once your simulation has produced sufficient events, the histograms written by the analysis manager form your simulated energy spectra. These spectra can then be visualized or further analyzed, for example to extract peak positions, resolutions, and efficiencies in later chapters that cover detailed analysis and detector effects.

Views: 8

Comments

Please login to add a comment.

Don't have an account? Register now!