将高斯拟合到 ROOT 直方图
Fitting gaussians to a ROOT histogram
我正在尝试编写一个程序来将多个高斯拟合到 ROOT 直方图,但不幸的是我对 pyROOT 的经验不足。
我有一个 Ba133 发射光谱的直方图,我想对峰进行高斯分布,以便我知道所述峰的 x 轴值,这是为了校准检测器。理想情况下,程序会迭代并找到峰值本身,但我可以自己指定它们。
目前我唯一的代码是:
import ROOT
infile = ROOT.TFile("run333.root")
Ba_spectrum = infile.Get("hQ0")
Ba_spectrum.Draw()
如果有人可以告诉我如何使用 pyroot 将高斯拟合到这些峰,最好是自动化的,将不胜感激。
谢谢
鉴于获得合适的拟合结果通常取决于从合理的初始参数值开始,您最好指定峰值的大致位置。例如,您可以有一个文本文件,其中包含对所有明显峰的高度、均值和宽度的猜测。
16000.0 625.0 25.0
500.0 750.0 50.0
...
然后运行像这样合身。
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))
我正在尝试编写一个程序来将多个高斯拟合到 ROOT 直方图,但不幸的是我对 pyROOT 的经验不足。
我有一个 Ba133 发射光谱的直方图,我想对峰进行高斯分布,以便我知道所述峰的 x 轴值,这是为了校准检测器。理想情况下,程序会迭代并找到峰值本身,但我可以自己指定它们。
目前我唯一的代码是:
import ROOT
infile = ROOT.TFile("run333.root")
Ba_spectrum = infile.Get("hQ0")
Ba_spectrum.Draw()
如果有人可以告诉我如何使用 pyroot 将高斯拟合到这些峰,最好是自动化的,将不胜感激。
谢谢
鉴于获得合适的拟合结果通常取决于从合理的初始参数值开始,您最好指定峰值的大致位置。例如,您可以有一个文本文件,其中包含对所有明显峰的高度、均值和宽度的猜测。
16000.0 625.0 25.0
500.0 750.0 50.0
...
然后运行像这样合身。
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))