Fitting gaussians to a ROOT histogram

Question

I'm trying to write a program to fit several gaussians to a ROOT histogram, but unfortunately my inexperience with pyROOT is showing.

I have a histogram of an emission spectrum of Ba133, and would like to fit gaussians to the peaks so that I know the x-axis value for said peaks, this in order to calibrate a detector. Ideally the program would iterate along and find the peaks itself, but I'm fine with having to specify them myself.

Currently the only code I have is:

import ROOT

infile = ROOT.TFile("run333.root")

Ba_spectrum = infile.Get("hQ0")

Ba_spectrum.Draw()

If someone could please tell me how to use pyroot to fit gaussians to these peaks, preferably automated, it would be greatly appreciated.

Thanks

Answer 1

Given that getting a decent fit result often depends on starting with reasonable initial parameter values, you're probably better off specifying the approximate locations of the peaks to begin with. You could for instance have a text file containing guesses of the heights, means, and widths for all the apparent peaks.

16000.0 625.0 25.0
  500.0 750.0 50.0
...

Then run the fits like this.

import ROOT

in_file = ROOT.TFile("run333.root")
if not in_file.IsOpen():
    raise SystemExit("Could not open file.")

spectrum = in_file.Get("hQ0")
if spectrum is None:
    raise SystemExit("Could not find spectrum histogram.")

N_GAUSS_PARAMS = 3

init = []
with open("init.txt") as f:
    for s in f:
        tokens = s.split()
        if not tokens:
            continue
        if len(tokens) != N_GAUSS_PARAMS:
            raise SystemExit(
                "Bad line in initial-value file: \"{}.\"".format(s)
            )

        init.append([float(t) for t in tokens])

n_peaks  = len(init)
n_params = N_GAUSS_PARAMS * n_peaks

fit_function = ROOT.TF1(
    "fit_function",
    "+".join(
        ["gaus({})".format(i)
         for i in range(0, n_params, N_GAUSS_PARAMS)]
    ), 0.0, 4100.0
)
for i in range(n_peaks):
    for j in range(N_GAUSS_PARAMS):
        fit_function.SetParameter(i * N_GAUSS_PARAMS + j, init[i][j])

spectrum.Fit(fit_function)

for i in range(1, n_params, N_GAUSS_PARAMS):
    print(fit_function.GetParameter(i))