Problems in fitting with ROOT 6.32

Roofit and RooStats
1

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

immagine
immagine779×537 31.9 KB

the red is the fit.
(With the previous ROOT version I had installed 6.30, everything worked just fine :smiling_face_with_tear:)

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

2

Hi Elia,

Thanks for the question.
I need to defer this to @jonas as you seem to point to a possible regression between 6.30 and 6.32.

Cheers,
D

3

Hi @Elia_Giulio_Grandoni, thanks for giving ROOT 6.32 a shot!

Let’s try to avoid this non-standard memory management first. Maybe your fit problem is also related to it.

The problem is that the lifetime of your RooDataHists, is not managed: they only live as long as the iteration in your loop, but the RooHistPdfs depend on them. can you also put them in a list so they survive? I think you don’t need the deep-copy then:

datahists_minv = []
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 = ROOT.RooHistPdf(f"pdf_minv_{name}", f"pdf_minv_{name}", ROOT.RooArgSet(minv_var), roohist_minv)

    datahists_minv.append(roohist_minv)
    pdf_minv.append(pdf_minv_one)

Let me know if this fixes the problem already!

Cheers,
Jonas

1 Like
4

Thanks, @jonas!

That magically solves the problem :tada:
Interestingly, this issue didn’t exist in the previous ROOT version.

Cheers,
Elia

5

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.