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

Would [TSpectrum](https://root.cern.ch/doc/master/classTSpectrum.html) work in your case? or you know the rough position of the peaks and you want to fit their parameters? — user2148414, Feb 09 '17 at 06:28
Potentially, like I said I've very little experience with ROOT, so I'm open to suggestions. I know the rough positions of the peaks simply by eyeballing them, I need to fit Gaussians to determine the x value, so that e.g. I can determine that a peak at 375 keV corresponds to bin no. 2500 for example. — ponsly, Feb 11 '17 at 15:43

Sam Marinelli · Accepted Answer · 2017-08-17T21:23:29.703

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))

That looks ideal, forgive my ignorance, but how do I then get the actual means out from this code? Is there a way to print or write them to a file? Thanks — ponsly, Feb 22 '17 at 15:18

Fitting gaussians to a ROOT histogram

1 Answers1