Failed Hypatia2 toy fits

[Dongliang]

Dear experts,
I’m using RooHypatia2 to fit a distribution and I get “HESSE STATUS=NOT POSDEF” and “MIGRAD STATUS=FAILED”. To investigate the issue I try some toy study using RooMCStudy and find that there are 5-15% of the fits will fail even if I generate 50k toy events. The script I used to test for 1k toy events is here:

#!/usr/bin/env python3
import ROOT

def check1(nEvt=1000, nTest=500, outfile = 'test_Jul11a.root'):
    w = ROOT.RooWorkspace('w')
    w.factory('Hypatia2::pdf(x[5610,5520,5720],l[-2.83,-7,-1],zeta[0,0,10],fb[0,-1,1],sigma[13.,0.1,30],mean[5620],aL[1.023,0.1,9],nL[1.565,0.1,9],aR[1.612,0.1,12],nR[1.507,0.1,12])')
    w['zeta'].setConstant(True)
    w['fb'].setConstant(True)
#     w['l'].setConstant(True)
#     w['sigma'].setConstant(True)
#     w.factory('CrystalBall::pdf(x[5610,5520,5720],mean[5620],sigma[13.,0.1,30],sigma,aL[1.023,0.1,9],nL[1.565,0.1,9],aR[1.612,0.1,12],nR[1.507,0.1,12])')
    obs = ROOT.RooArgSet(w['x'])
    pdf = w.pdf('pdf')

    mcstudy = ROOT.RooMCStudy(pdf, obs, ROOT.RooFit.FitOptions(ROOT.RooFit.Save(True)))
    mcstudy.generateAndFit(nTest,nEvt);

    nGood = 0
    fout1 = ROOT.TFile(outfile,"recreate")
    fout1.cd()
    for i in range(nTest):
#         mcstudy.fitResult(i).Print()
        if mcstudy.fitResult(i).status()==0: nGood += 1
        mcstudy.fitResult(i).Write(mcstudy.fitResult(i).GetName()+f"_{i}")
    fout1.Close()

    print(f"{nGood}/{nTest} fits are good.")
    a = input('wait...')
check1()

The fraction of successful fits is printed out at the end. I get 440/500 in my test. Is this normal? Thanks.

PS, I also tried RooCrystalBall for comparison, the toy fits of which all converge without error. In case it’s useful, I also have another function to check the parameters from those failed fits
check_toy_study1.py (5.2 KB). It seems that in most of the failed fits, the parameters (red) are similar to the successful fits (blue). And I’m using LCG_102/ROOT/6.26.04/x86_64-centos7-gcc11-opt.

[bellenot]

@jonas can most probably help you

[jonas]

Hi @Dongliang,

again a Hypathia2 question! I hope I can answer this one better this time I hope I can also do a detailed investigation on the normalization problem that you also reported on the forum.

The reason why the crystal ball converges better than the RooHypathia2 is that the Hypathia2 doesn’t support analytic integration yet (although there is a GitHub issue open about this that I also still need to work on: [RF] RooHypatia2 Analytical integral integration · Issue #7254 · root-project/root · GitHub).

Until then, the best you can do is to play around with the numeric integrator and also some Minuit options to increase the probability for convergence. For example, I get 500/500 good fits with Minuit strategy setting “2” and the “RooAdaptiveGaussKronrodIntegrator1D” with stricter tolerance:

#!/usr/bin/env python3
import ROOT

ROOT.RooAbsReal.defaultIntegratorConfig().setEpsAbs(1e-9)
ROOT.RooAbsReal.defaultIntegratorConfig().setEpsRel(1e-9)
ROOT.RooAbsReal.defaultIntegratorConfig().method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")

def check1(nEvt=1000, nTest=500, outfile = 'test_Jul11a.root'):
    w = ROOT.RooWorkspace('w')
    w.factory('Hypatia2::pdf(x[5610,5520,5720],l[-2.83,-7,-1],zeta[0,0,10],fb[0,-1,1],sigma[13.,0.1,30],mean[5620],aL[1.023,0.1,9],nL[1.565,0.1,9],aR[1.612,0.1,12],nR[1.507,0.1,12])')
    w['zeta'].setConstant(True)
    w['fb'].setConstant(True)
#     w['l'].setConstant(True)
#     w['sigma'].setConstant(True)
#     w.factory('CrystalBall::pdf(x[5610,5520,5720],mean[5620],sigma[13.,0.1,30],sigma,aL[1.023,0.1,9],nL[1.565,0.1,9],aR[1.612,0.1,12],nR[1.507,0.1,12])')
    obs = ROOT.RooArgSet(w['x'])
    pdf = w.pdf('pdf')

    mcstudy = ROOT.RooMCStudy(pdf, obs, FitOptions=dict(Save=True, Strategy=2))
    mcstudy.generateAndFit(nTest,nEvt);

    nGood = 0
    fout1 = ROOT.TFile(outfile,"recreate")
    fout1.cd()
    for i in range(nTest):
#         mcstudy.fitResult(i).Print()
        if mcstudy.fitResult(i).status()==0: nGood += 1
        mcstudy.fitResult(i).Write(mcstudy.fitResult(i).GetName()+f"_{i}")
    fout1.Close()

    print(f"{nGood}/{nTest} fits are good.")
    # a = input('wait...')

check1()

Here is a tutorial that tells you more about configuring the numeric integrator: ROOT: tutorials/roofit/rf901_numintconfig.py File Reference.

I hope this helps you with your fit!

Jonas

[Dongliang]

Hi @jonas , thanks a lot for the reply. This solves my problem perfectly!
I came across the GitHub issue, and I think the analytical integration could be enabled with minor modification to the current (commented-out) code. We can continue the discussion there…