23.10. Analyzing the Spectrum with ROOT

Table of Contents

Preparing the ROOT Environment

To analyze the gamma ray detector spectrum, you will use the ROOT framework to read and plot the data that your Geant4 application produced. At this stage, you should already have a ROOT file created by G4AnalysisManager that contains at least one 1D histogram of deposited energy in the scintillator.

Start by launching ROOT from a terminal using the root command. You can work either in the interactive C++ prompt (the ROOT shell, called CINT or Cling in recent ROOT versions) or by writing short macros in .C files. For beginners, the interactive prompt is usually more convenient.

Once ROOT is open, you will load your output file with the TFile class. ROOT works with objects stored in files, and you can list the contents, retrieve histograms, and immediately draw them on a canvas for inspection.

It is good practice to know the exact names of your histograms and ntuples from your Geant4 code. The names you pass to CreateH1, CreateNtuple, and similar methods in G4AnalysisManager will appear as object names in the ROOT file.

Always make sure that the ROOT file you want to analyze is closed properly by Geant4 before opening it in ROOT. Opening a file that is still being written can lead to corrupted or incomplete data.

Opening Geant4 Output in ROOT

Assume that in your gamma ray detector example, you configured the analysis manager to write a file such as gammaDetector.root and created a 1D histogram of energy deposition called "Edep" in a directory called "h1".

To open the file in ROOT, start ROOT and run:

cpp

TFile *f = TFile::Open("gammaDetector.root");

You can check what is inside by listing the top level content:

cpp

f->ls();

If you used the default G4AnalysisManager convention with one dimensional histograms, you may see an object directory named "h1" that contains your histograms. You can navigate into it:

cpp

f->cd("h1");
gDirectory->ls();

Now you should see your histogram, for example "Edep". Retrieve it as a TH1 object:

cpp

TH1D *hE = (TH1D*)gDirectory->Get("Edep");

If you are not sure about the type, you can use TH1 *:

cpp

TH1 *hE = (TH1*)gDirectory->Get("Edep");

At this point, hE contains your energy spectrum and is ready to be drawn, rebinned, fitted, or otherwise processed.

If your Geant4 application wrote ntuples instead of or in addition to histograms, you will also see TTree objects in the ROOT file, often named "ntuple" or similar. You can list them with:

cpp

f->ls();

and then:

cpp

TTree *t = (TTree*)f->Get("ntuple");
t->Print();

This prints the branches and their types. For a gamma ray detector, you might have branches like energy deposition, detector ID, or event ID.

When you use Get to retrieve objects from a ROOT file, always check that the returned pointer is not nullptr before using it. A null pointer indicates that the object name or path is wrong.

Drawing and Inspecting the Spectrum

Once you have the histogram in memory, you can inspect the gamma ray energy spectrum visually. ROOT uses a TCanvas to display graphics.

Create a new canvas and draw the histogram:

cpp

TCanvas *c1 = new TCanvas("c1", "Gamma-ray energy spectrum", 800, 600);
hE->Draw();

By default, ROOT draws a step-like histogram. You can change the draw option, for example:

cpp

hE->Draw("HIST");
hE->SetLineColor(kBlue);
hE->SetLineWidth(2);

You can also change axis titles so that the plot is physically meaningful:

cpp

hE->GetXaxis()->SetTitle("Deposited energy [MeV]");
hE->GetYaxis()->SetTitle("Counts");
c1->Update();

Zooming into a region of interest can be done interactively with the mouse, or programmatically by setting axis ranges:

cpp

hE->GetXaxis()->SetRangeUser(0.0, 2.0);  // example: 0 to 2 MeV
c1->Modified();
c1->Update();

If your Geant4 histogram was created with units in mind (for example, filling with edep / MeV), then the x axis is already in MeV. Double check the filling code so that your axis labels correspond to the correct units.

You may find that the statistics box in the top right corner shows total entries, mean, and RMS of the energy distribution. To make it visible or adjust options:

cpp

gStyle->SetOptStat(1110);  // show entries, mean, RMS
c1->Modified();
c1->Update();

Always confirm that the units used to fill the histogram in Geant4 match the units you display in ROOT. A common mistake is to fill energies in keV but label the axis as MeV, which leads to misinterpretation of spectral features.

Basic Statistical Quantities

ROOT can provide basic statistical information about your gamma ray spectrum through methods of TH1. For instance, you can quickly inspect global properties:

cpp

Int_t entries = hE->GetEntries();
Double_t mean   = hE->GetMean();
Double_t rms    = hE->GetRMS();

GetMean and GetRMS operate in the x variable, so here they represent the mean deposited energy and the width of the distribution. These global values are useful for a quick check, but for gamma spectroscopy you are usually more interested in specific peaks rather than the entire distribution.

You can also compute integrals over selected ranges in order to count events in particular energy windows. If you know the bin numbers that correspond to a range, use:

cpp

Int_t bin_low  = hE->FindBin(0.45);  // example: lower energy in MeV
Int_t bin_high = hE->FindBin(0.55);  // example: upper energy in MeV
Double_t counts = hE->Integral(bin_low, bin_high);

Integral returns the number of entries in that bin range without any corrections. You can divide this by the total number of simulated primary particles to estimate detection probabilities or efficiencies.

When you work with histograms that have variable bin widths, be aware that Integral can be normalized by bin width. For a uniform binning in this simple example, a plain integral is enough.

When interpreting counts in a given peak or energy region, always normalize by the number of primary gammas simulated if you want to compare simulations that have different event counts or if you compare with experimental efficiencies.

Identifying and Fitting Photopeaks

A central goal of analyzing a gamma ray detector spectrum is to identify characteristic features such as photopeaks, Compton edges, and backscatter peaks. In a simple example with a monoenergetic gamma source, the main visible structure is typically the photopeak at or near the input gamma energy.

To identify the photopeak, visually inspect the spectrum in ROOT and locate the bin corresponding to the maximum of the peak:

cpp

Int_t maxBin = hE->GetMaximumBin();
Double_t E_peak_estimate = hE->GetBinCenter(maxBin);

This gives an approximate peak position. To obtain a more precise value and an estimate of the detector energy resolution, you can fit the peak with a Gaussian function. In ROOT, use TF1 to define a Gaussian and fit it to a selected energy range around the peak.

For example, suppose you expect a photopeak around 0.662 MeV (from a Cs-137 source). You can define a fit range around that energy:

cpp

TF1 *gaus = new TF1("gaus", "gaus", 0.5, 0.8);  // fit between 0.5 and 0.8 MeV
hE->Fit(gaus, "R");  // "R" means use the defined range

The Gaussian function in ROOT has three parameters. After the fit, you can access them:

cpp

Double_t norm   = gaus->GetParameter(0);  // amplitude
Double_t mean   = gaus->GetParameter(1);  // peak position
Double_t sigma  = gaus->GetParameter(2);  // standard deviation
Double_t meanErr  = gaus->GetParError(1);
Double_t sigmaErr = gaus->GetParError(2);

The standard deviation is related to the full width at half maximum by:
$$ \text{FWHM} = 2.355 \, \sigma $$
You can compute it in ROOT:

cpp

Double_t fwhm = 2.355 * sigma;
Double_t resolution = fwhm / mean;   // relative energy resolution

Remember the relation between Gaussian width and resolution:
$$ \boxed{\text{FWHM} = 2.355 \, \sigma} $$
and energy resolution at energy $E$ is often expressed as:
$$ \boxed{R = \frac{\text{FWHM}}{E}} $$
where $R$ is the relative resolution.

You can print or store these values for later comparison with experimental data or with other simulations. This is the main bridge from raw simulation output to detector performance metrics.

If your spectrum shows more than one peak, you can restrict the fit range to each peak individually, or use more complex composite functions that sum several Gaussians plus a background term.

Comparing Ideal and Smeared Spectra

In a simple Geant4 gamma detector example, you might first record an "ideal" energy spectrum where the deposited energy is filled without any detector resolution effects. In a later step, you may apply energy smearing (for example, Gaussian smearing) either in Geant4 or in ROOT to simulate a realistic detector response.

To compare an ideal spectrum and a smeared one, open the ROOT file that contains both histograms. For example, you might have "Edep_ideal" and "Edep_smeared".

Retrieve both:

cpp

TH1 *hIdeal   = (TH1*)gDirectory->Get("Edep_ideal");
TH1 *hSmeared = (TH1*)gDirectory->Get("Edep_smeared");

Normalize them to the same area if they have different total entries, so that the shapes can be compared directly:

cpp

hIdeal->Scale(1.0 / hIdeal->Integral());
hSmeared->Scale(1.0 / hSmeared->Integral());

Now draw them on the same canvas:

cpp

TCanvas *c2 = new TCanvas("c2", "Ideal vs smeared spectra", 800, 600);
hIdeal->SetLineColor(kBlack);
hIdeal->SetLineWidth(2);
hIdeal->Draw("HIST");
hSmeared->SetLineColor(kRed);
hSmeared->SetLineWidth(2);
hSmeared->Draw("HIST SAME");

You should observe that in the smeared spectrum, photopeaks become broader and possibly lower in amplitude while the total area remains the same after normalization. This visual comparison helps you understand how detector resolution affects peak shapes and the ability to resolve nearby energy lines.

If you did not implement smearing inside Geant4, you can apply it in ROOT by creating a new histogram and filling it with smeared values. Loop over the original histogram bins and for each bin center and content, sample a Gaussian with given $\sigma$ and fill the new histogram. For many beginners, however, implementing smearing on the Geant4 side is simpler and more consistent with later analysis.

When comparing ideal and smeared spectra, always normalize the histograms before visually comparing shapes. Comparing raw counts can be misleading when the total number of events or binning differs between spectra.

Saving Plots and Results

Once you have an informative spectrum plot and fitted parameters, you will usually want to save the plot and possibly some numerical results.

To save the canvas as an image, use:

cpp

c1->SaveAs("gamma_spectrum.png");    // or ".pdf", ".eps", ".root"

You can do the same for comparison plots:

cpp

c2->SaveAs("ideal_vs_smeared.pdf");

If you want to preserve the ROOT objects for later, you can write them into a new ROOT file:

cpp

TFile *fout = TFile::Open("analysisResults.root", "RECREATE");
c1->Write("c_spectrum");
c2->Write("c_compare");
gaus->Write("fit_photopeak");
fout->Close();

You can also save numerical values such as fitted peak position and resolution into a text file directly from ROOT:

cpp

std::ofstream txt("peak_results.txt");
txt << "Peak mean [MeV]: " << mean << "\n";
txt << "FWHM [MeV]: " << fwhm << "\n";
txt << "Resolution (FWHM/E): " << resolution << "\n";
txt.close();

This gives you a compact record of your main analysis results that you can later use for reports or further comparison, for example when you change detector geometry or materials in the Geant4 simulation.

Always keep a clear link between the ROOT analysis and the Geant4 configuration that produced the data. Record the simulation version, geometry, materials, physics list, and source settings together with your analysis results to ensure that your findings can be reproduced.

By following these steps, you turn the raw energy deposition data from your Geant4 gamma ray detector simulation into interpretable spectra and quantitative measures of detector performance using ROOT.