Dear experts,
I recently downloaded the new root version and my code is doing strange things in the fitting. It is a simple model signal + background and I’m trying to fit the number of signal events from a dataset generated from the RooHistPdf from imported histograms.
The code is here
import ROOT
import numpy as np
import globalConstants as ctes
from copy import deepcopy as cp
# Variables
minv_var = ROOT.RooRealVar("m2inv", "m^{2}_{inv}", -200, 500, "[MeV^{2}]")
minv_var.setBins(ctes.BINS_OBS_MASS)
pdf_minv = []
for i in range(3):
name = ctes.NAMES[i]
# Input file
inputF = ROOT.TFile(f"{ctes.FILE_PATH}/Histo{name}Fit.root")
# Histograms from input file
histo_minv = inputF.Get("himas")
histo_minv.Rebin(ctes.BINNING_HISTO_MASS)
# RooHist to PDF
roohist_minv = ROOT.RooDataHist(f"histo_minv_{name}", f"histo_minv_{name}", ROOT.RooArgList(minv_var), histo_minv)
pdf_minv_one = cp(ROOT.RooHistPdf(f"pdf_minv_{name}", f"pdf_minv_{name}", ROOT.RooArgSet(minv_var), roohist_minv))
pdf_minv.append(pdf_minv_one)
'''
c = ROOT.TCanvas()
frame = minv_var.frame()
pdf_minv_one.plotOn(frame)
frame.Draw()
c.SaveAs(f"{ctes.IMMAGINI_PATH}/pdf{name}.pdf")
'''
nevBKG = ctes.NRMD + ctes.NACC
nevSIG = np.round(nevBKG * ctes.BR)
evTOT = (nevBKG + nevSIG)
NRMD_var = ROOT.RooRealVar("NRMD_var", "NRMD_var", ctes.NRMD, 0.75*ctes.NRMD, 1.25*ctes.NRMD)
NACC_var = ROOT.RooRealVar("NACC_var", "NACC_var", ctes.NACC, 0.75*ctes.NACC, 1.25*ctes.NACC)
NSIG_var = ROOT.RooRealVar("NSIG_var", "NSIG_var", nevSIG, -4000, 4000)
'''
# Create model
model_mass = ROOT.RooAddPdf("model_mass", "model_mass",
ROOT.RooArgList(pdf_minv[0], pdf_minv[1], pdf_minv[2]),
ROOT.RooArgList(NRMD_var, NACC_var, NSIG_var))
'''
model_mass = ROOT.RooAddPdf("model_mass", "model_mass",
ROOT.RooArgList(pdf_minv[0], pdf_minv[1], pdf_minv[2]),
ROOT.RooArgList(NRMD_var, NACC_var, NSIG_var))
# Generate data
data_minv = model_mass.generateBinned(ROOT.RooArgSet(minv_var), evTOT)
data_minv.SetName("data_minv")
res = model_mass.fitTo(data_minv, ROOT.RooFit.Save())
res.Print("v")
c = ROOT.TCanvas()
frame_minv = minv_var.frame()
data_minv.plotOn(frame_minv)
model_mass.plotOn(frame_minv, ROOT.RooFit.LineColor(ROOT.kRed))
model_mass.plotOn(frame_minv, ROOT.RooFit.Components(pdf_minv[0]))
model_mass.plotOn(frame_minv, ROOT.RooFit.Components(pdf_minv[1]))
model_mass.plotOn(frame_minv, ROOT.RooFit.Components(pdf_minv[2]))
frame_minv.Draw()
c.SaveAs(f"{ctes.IMMAGINI_PATH}minv_fit.pdf")
The fit has a minimum for -17k signal events…that does not make any sense. If a fix the lower bound at higher values it just converges for that value (but it still does not make any sense). With a lower bound of -3k signal events this is the output of the fit
[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Fitting -- RooAbsPdf::fitTo(model_mass) fixing normalization set for coefficient determination to observables in data
[#1] INFO:Fitting -- using CPU computation library compiled with -mavx512
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_model_mass_data_minv) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
Minuit2Minimizer: Minimize with max-calls 1500 convergence for edm < 1 strategy 1
Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator
Info in <Minuit2>: MnSeedGenerator Initial state: FCN = -36170.34841 Edm = 1658833.674 NCalls = 13
Info in <Minuit2>: MnSeedGenerator Initial state
Minimum value : -36170.34841
Edm : 1658833.674
Internal parameters: [ 0 0 0]
Internal gradient : [ 793.9256841 -794.5498004 2983.818195]
Internal covariance matrix:
[[ 0.54667118 0 0]
[ 0 0.0018598228 0]
[ 0 0 0.70644312]]]
Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 1500
Info in <Minuit2>: VariableMetricBuilder 0 - FCN = -36170.34841 Edm = 1658833.674 NCalls = 13
Info in <Minuit2>: VariableMetricBuilder 1 - FCN = -38717.22217 Edm = 965123.3908 NCalls = 29
Info in <Minuit2>: VariableMetricBuilder 2 - FCN = -39006.62401 Edm = 68632.75995 NCalls = 45
Info in <Minuit2>: VariableMetricBuilder 3 - FCN = -39599.43306 Edm = 44136.29004 NCalls = 53
Info in <Minuit2>: VariableMetricBuilder 4 - FCN = -40149.78342 Edm = 131.8522205 NCalls = 61
Info in <Minuit2>: VariableMetricBuilder 5 - FCN = -40263.46092 Edm = 6.186113627 NCalls = 68
Info in <Minuit2>: VariableMetricBuilder 6 - FCN = -40275.4475 Edm = 1.527982278 NCalls = 77
Info in <Minuit2>: VariableMetricBuilder 7 - FCN = -40277.76785 Edm = 0.1074850123 NCalls = 85
Info in <Minuit2>: VariableMetricBuilder 8 - FCN = -40277.96859 Edm = 0.001630252139 NCalls = 93
Info in <Minuit2>: VariableMetricBuilder 9 - FCN = -40277.97102 Edm = 7.042035345e-06 NCalls = 101
Info in <Minuit2>: VariableMetricBuilder After Hessian
Info in <Minuit2>: VariableMetricBuilder 10 - FCN = -40277.97102 Edm = 2.725351439e-06 NCalls = 119
Minuit2Minimizer : Valid minimum - status = 0
FVAL = -40277.9710240491622
Edm = 2.7253514392615243e-06
Nfcn = 119
NACC_var = 2428.5 +/- 0.517524 (limited)
NRMD_var = 17208.4 +/- 130.98 (limited)
NSIG_var = -3000 +/- 0.504051 (limited)
Info in <Minuit2>: Minuit2Minimizer::Hesse Using max-calls 1500
Info in <Minuit2>: Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
RooFitResult: minimized FCN value: -40278, estimated distance to minimum: 2.71947e-06
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter InitialValue FinalValue +/- Error GblCorr.
-------------------- ------------ -------------------------- --------
NACC_var 3.2380e+03 2.4285e+03 +/- 5.18e-01 <none>
NRMD_var 1.4029e+04 1.7208e+04 +/- 1.31e+02 <none>
NSIG_var 0.0000e+00 -3.0000e+03 +/- 5.04e-01 <none>
It should converge around zero since the dataset has zero signal events…
Here is an image of the fit result
the red is the fit.
(With the previous ROOT version I had installed 6.30, everything worked just fine )
Another unrelated question is: Why if I do not use
deepcopy
for the pdf the program crashes? It seems there are some problems with memory handling…
Thanks in advance for the advices,
Elia
_ROOT Version:6.32.00