Hi everyone,
I’m a beginner trying to learn more about ROOT and data analysis. From a course on radioactivity, I learned the basic principles of gamma spectroscopy (see attached photo), and for fun I decided to try calibrating the energy scale of a detector using recognizable photopeaks.
To do this, I need to identify the peak positions by fitting them with a model that includes a Gaussian (for the peak) on top of a linear background.
However, I’m having trouble performing a good fit using RooFit. Specifically, I’m trying to model the data with a linear + Gaussian PDF, but the fit I obtain looks very poor compared to the data (see attached plot). On the other hand, using a simple TF1 fit gives a result that matches the data much better.
I attach the code I used. I’ve tried two approaches:
- Fitting the full model (linear + Gaussian) at once, with all parameters floating.
- Fitting just the background (intercept and slope in the sidebands), fixing them afterwards, and then fitting the Gaussian on top with only the signal and background normalizations (
nsig,nbkg) floating.
In both cases, the result looks far from optimal, and I’m not sure what I’m doing wrong.
I’d really appreciate any advice or suggestions!
Thanks in advance,
Luca
import numpy as np
import ROOT
from ROOT import RooFit, RooRealVar, RooDataHist, RooArgList, RooGenericPdf, RooGaussian, RooAddPdf, RooFormulaVar
data21 = np.loadtxt("M21_00000list.txt", comments="#", skiprows=3)
channels = data21[:, 1]
h_data21 = ROOT.TH1F("h_data21", "Spettro M21", 2048, 0, 2048)
for c in channels:
h_data21.Fill(c)
channel = RooRealVar("channel", "channel", 0, 2048)
rh_data21 = RooDataHist('rh_data21','rh_data21',RooArgList(channel),h_data21)
channel.setRange('left',700,740)
channel.setRange('right',860,900)
channel.setRange('all',700,900)
slope = RooRealVar('slope','slope',-.1,-10,0.)
offset = RooRealVar('offset','offset',1000,0,2000)
linear = RooGenericPdf('linear','offset + slope * channel',RooArgList(offset,slope,channel))
#linear.fitTo(rh_data21, RooFit.Range('left,right'))
#slope.setConstant(1)
#offset.setConstant(1)
mean = RooRealVar('mean','mean', 810, 800, 825)
sigma = RooRealVar('sigma','sigma',10,0,40)
gauss = RooGaussian('gauss','gauss',channel,mean,sigma)
nsig = RooRealVar('nsig','nsig',1000,100,8000)
nbkg = RooRealVar('nbkg','nbkg',10000,100,30000)
#nbkg = RooFormulaVar('nbkg','nbkg','-20000*slope',RooArgList(slope))
#nbkg.setConstant(1)
pdftot = RooAddPdf('pdftot','pdftot',RooArgList(gauss,linear),RooArgList(nsig,nbkg))
pdftot.fitTo(rh_data21, RooFit.Range('all'))
frame21 = channel.frame()
rh_data21.plotOn(frame21,RooFit.MarkerStyle(7))
#linear.plotOn(frame21, RooFit.LineColor(ROOT.kRed), RooFit.LineWidth(2))
pdftot.plotOn(frame21, RooFit.LineColor(ROOT.kBlue), RooFit.LineWidth(3))
c21 = ROOT.TCanvas('c21','c21')
c21.SetLogy()
frame21.GetXaxis().SetRangeUser(600,1000)
frame21.Draw()
c21.SaveAs('c21.svg')
(total spectrum)
(ROOT TF1 Fit)
(RooFit)
