35.10. Analyze the Results with ROOT

Table of Contents

Preparing ROOT for Final Project Analysis

For the final project, you already have a Geant4 simulation that writes output in a ROOT compatible format using G4AnalysisManager. This chapter focuses only on what you do after the simulation has finished: how to open the ROOT file, inspect its contents, make plots, apply cuts, and extract physics results.

You will work mainly inside a ROOT session or a small ROOT macro, not inside Geant4.

Inspecting the Output ROOT File

Once your Geant4 run is complete, you should have one or more ROOT files, for example finalProject.root. Start an interactive ROOT session from the command line:

bash

root finalProject.root

ROOT automatically opens a prompt where you can inspect the file content.

Use TFile::ls() to see what objects are stored:

cpp

TFile *f = (TFile*)gROOT->GetListOfFiles()->At(0);
f->ls();

If you used G4AnalysisManager to create histograms and ntuples, you will typically see:

One or more TH1 or TH2 objects, such as hEdep.
One TTree or TNtuple, such as ntuple or Ntuple0.

To explore the branches of a TTree, for example:

cpp

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

This prints branch names and types, such as Edep, x, y, z, time, eventID, detID and any other columns you defined.

Always verify branch names and types in the ROOT file before writing analysis code. Using a wrong branch name or assuming a wrong type often leads to silent mistakes in your results.

If you used G4AnalysisManager’s default naming, ntuples may be called "ntuple0", "ntuple1", and histograms "h0", "h1", and so on. Check your Geant4 code or the ROOT file listing to find the exact identifiers.

Plotting Basic Distributions

ROOT can quickly plot histograms and tree variables. If your Geant4 application created a histogram, such as an energy spectrum called Edep, you can draw it directly:

cpp

TH1D *hEdep = (TH1D*)f->Get("Edep");
hEdep->Draw();

If you only stored an ntuple (tree) with columns, use the ROOT Draw method:

cpp

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

ROOT automatically creates a histogram of the variable Edep using default binning. To specify the number of bins and range:

cpp

t->Draw("Edep>>hE(200, 0, 2.0)"); // 200 bins from 0 to 2.0 (e.g. MeV)
TH1D *hE = (TH1D*)gDirectory->Get("hE");

You can set axis titles and style for clearer plots:

cpp

hE->GetXaxis()->SetTitle("Deposited energy [MeV]");
hE->GetYaxis()->SetTitle("Counts");
gStyle->SetOptStat(1110); // Show entries, mean, RMS
hE->Draw();

Use logy for spectra that span several orders of magnitude:

cpp

gPad->SetLogy();
hE->Draw();

You can quickly check other distributions in the same way, for example position or time:

cpp

t->Draw("x");                 // X position distribution
t->Draw("time");              // Time distribution
t->Draw("y:x>>hXY(100, -50, 50, 100, -50, 50)"); // 2D hit map

For 2D histograms, choose a color representation:

cpp

TH2D *hXY = (TH2D*)gDirectory->Get("hXY");
hXY->Draw("COLZ");

These quick plots help you visually confirm that your Geant4 output is reasonable before doing more detailed analysis.

Always start with simple sanity checks: total entries, typical energy scale, and basic shapes. This often reveals unit mistakes or incorrect cuts very early.

Applying Cuts and Selections

To extract physics quantities such as detector efficiency or background rejection, you usually need to apply selections on the recorded variables. ROOT lets you express cuts as simple C-like strings.

For example, to plot energy only in the central detector, assuming a branch detID:

cpp

t->Draw("Edep>>hCenter(200,0,2.0)", "detID == 0");

To select events above a threshold, for example only hits with Edep > 0.1 MeV:

cpp

t->Draw("Edep>>hAbove(200,0,2.0)", "Edep > 0.1");

To combine cuts:

cpp

t->Draw("Edep>>hSel(200,0,2.0)",
        "detID == 0 && Edep > 0.1 && time < 50");

For coincidence-like conditions that depend on multiple hits in the same event, you often need a more detailed analysis loop, but you can still use ROOT selections to inspect subsets first.

A very useful feature is to compute statistics with cuts using TTree::Draw and TTree::GetSelectedRows. For example, to count the number of hits above a threshold:

cpp

Long64_t nAll  = t->GetEntries();
Long64_t nPass = t->Draw("Edep", "Edep > 0.1", "goff");
double efficiency = double(nPass) / double(nAll);

The "goff" option tells ROOT not to draw a histogram, only to use the selection internally.

When computing quantities like efficiency or transmission, always use the same definition consistently: clearly separate total incident events, detected events, and any energy or timing thresholds used in the selection.

Extracting Physics Quantities

The main goal of your final project analysis is to turn raw simulation output into physical results that answer your project’s scientific question. Typical quantities include energy spectra, efficiencies, depth dose curves, or spatial distributions.

Energy spectra and peak analysis

To analyze an energy spectrum, first create a clean histogram, then extract basic statistics and, if needed, fit peaks.

Assume you have a histogram hE with the total deposited energy per event in a detector. You can retrieve useful quantities:

cpp

int    entries = hE->GetEntries();
double mean    = hE->GetMean();
double rms     = hE->GetRMS();

To focus on a particular energy region, for example near a gamma peak, you can restrict the fit range and use a Gaussian model:

cpp

TF1 *gaus = new TF1("gaus", "gaus", 0.4, 0.6); // fit 0.4–0.6 MeV region
hE->Fit(gaus, "R");

After the fit, you can access fitted parameters:

cpp

double peakE   = gaus->GetParameter(1); // mean
double sigmaE  = gaus->GetParameter(2); // sigma
double fwhm    = 2.355 * sigmaE;
double res     = fwhm / peakE;          // relative energy resolution

You can then report energy resolution, compare to expectations, or compare between different detector configurations.

Efficiencies and probabilities

If you stored one entry per event with a total deposited energy EdepTotal, you can estimate the detection efficiency above some threshold:

cpp

TTree *tEvent = (TTree*)f->Get("eventNtuple");
Long64_t nTotal = tEvent->GetEntries();
Long64_t nDet   = tEvent->Draw("EdepTotal",
                               "EdepTotal > 0.2", "goff");
double eff      = double(nDet) / double(nTotal);
double sigmaEff = std::sqrt(eff * (1.0 - eff) / nTotal);

This efficiency with statistical uncertainty is usually a key result when comparing detector designs or materials.

Depth dose or spatial distributions

If your project involves a phantom or detector divided into bins (for example voxel ID or slice index), you probably stored an index and an energy per hit. To obtain a depth dose curve, you can accumulate energy as a function of bin index with a TProfile or histogram.

For example, suppose you have branches sliceID and Edep for each step or hit:

cpp

TTree *tHit = (TTree*)f->Get("hits");
TProfile *pDose = new TProfile("pDose","Depth dose",
                               100, 0, 100, 0, 1e6); // slice 0–99
tHit->Draw("Edep:sliceID>>pDose", "", "prof");
pDose->GetXaxis()->SetTitle("Slice index");
pDose->GetYaxis()->SetTitle("Mean deposited energy [MeV]");
pDose->Draw();

You can convert slice index to depth in mm or cm by setting custom axis labels or by using a histogram with an appropriate binning that reflects real distances.

Coincidences and correlations

If your final project involves coincidence detection or correlations between two detectors, you may have ntuples with one entry per hit and fields like eventID, detID, Edep, and time. In this case, a straightforward analysis is often easier with a small ROOT macro that loops over the tree, groups hits by event, and applies coincidence criteria.

For example, a simplified outline in a macro:

cpp

void analyzeCoincidences() {
  TFile *f = TFile::Open("finalProject.root");
  TTree *t = (TTree*)f->Get("hits");
  int eventID, detID;
  double Edep, time;
  t->SetBranchAddress("eventID", &eventID);
  t->SetBranchAddress("detID",   &detID);
  t->SetBranchAddress("Edep",    &Edep);
  t->SetBranchAddress("time",    &time);
  const double Ecut = 0.3;
  const double dtMax = 5.0; // time window, for example ns
  Long64_t nEntries = t->GetEntries();
  int currentEvent = -1;
  std::vector<int>    dets;
  std::vector<double> times;
  std::vector<double> edeps;
  int nCoinc = 0;
  for (Long64_t i = 0; i < nEntries; ++i) {
    t->GetEntry(i);
    if (eventID != currentEvent && currentEvent != -1) {
      // Analyze hits of the previous event
      for (size_t a = 0; a < dets.size(); ++a) {
        if (edeps[a] < Ecut) continue;
        for (size_t b = a+1; b < dets.size(); ++b) {
          if (edeps[b] < Ecut) continue;
          if (std::fabs(times[a] - times[b]) < dtMax) {
            ++nCoinc;
          }
        }
      }
      dets.clear();
      times.clear();
      edeps.clear();
    }
    currentEvent = eventID;
    dets.push_back(detID);
    times.push_back(time);
    edeps.push_back(Edep);
  }
  // Analyze last event
  // (repeat coincidence check here if needed)
  std::cout << "Total coincidences: " << nCoinc << std::endl;
}

You can adapt this structure to your final project needs, such as storing coincidence energy sums, line of response information, or other derived quantities.

When performing event-based analyses such as coincidences or track reconstruction, always make sure your tree is properly sorted by eventID. If not, group hits using both eventID and entry index, or sort them explicitly inside your macro.

Producing Publication-Quality Plots and Tables

For the final project report, you will likely need high quality figures and numerical summaries of your results. ROOT allows you to customize plots and export them to vector or bitmap formats.

To improve the appearance of a canvas:

cpp

TCanvas *c1 = new TCanvas("c1","Energy spectrum",800,600);
c1->SetGrid();
hE->SetLineColor(kBlue+1);
hE->SetLineWidth(2);
hE->SetTitle("Detector energy spectrum;Energy [MeV];Counts");
hE->Draw();
c1->SaveAs("energy_spectrum.pdf");
c1->SaveAs("energy_spectrum.png");

You can also combine multiple histograms in one plot to compare configurations, for example different materials or thicknesses:

cpp

TH1D *hLead  = (TH1D*)f->Get("Edep_lead");
TH1D *hAl    = (TH1D*)f->Get("Edep_aluminum");
hLead->SetLineColor(kRed);
hAl->SetLineColor(kBlue);
hLead->SetTitle("Energy spectra;Energy [MeV];Counts");
hLead->Draw();
hAl->Draw("SAME");
TLegend *leg = new TLegend(0.6,0.7,0.88,0.88);
leg->AddEntry(hLead, "Lead", "l");
leg->AddEntry(hAl,   "Aluminum", "l");
leg->Draw();
c1->SaveAs("comparison_materials.pdf");

To summarize key numerical results, you can print them directly from ROOT:

cpp

std::cout << "Detection efficiency: " << eff*100.0
          << " % ± " << sigmaEff*100.0 << " %" << std::endl;

You can copy these values into your report, or export them to a CSV file by writing them out yourself or using ROOT’s TTree::Scan or TTree::Draw("var >> file.csv") features.

Use consistent units and labels in all plots and tables, and explicitly state thresholds, cuts, and normalization choices. This is essential for others to interpret your final project results correctly.

Connecting Simulation and Analysis

At this stage of the course you should see a full chain:

Geant4 simulation produces ROOT data (ntuples, histograms).
ROOT is used to visualize the raw distributions.
Cuts and selections are applied to extract physically meaningful observables.
Fits, integrals, and counting give quantitative results with uncertainties.
Plots and tables summarize these results in a clear, reproducible way.

If you find during ROOT analysis that you are missing an important variable, or that your ntuple structure makes some analysis difficult, note this for later improvement of your Geant4 application. A well designed output format simplifies ROOT analysis, especially in larger projects.

For the final project, your goal is to demonstrate that you can complete this whole cycle: use simulation output, analyze it in ROOT, and present results that directly answer your chosen scientific question.