Creating own p.d.f. in RooFit

Roofit and RooStats
,
1

Dear experts,

I’m trying to create my own roofit p.d.f. with PyROOT to fit some distribution. It should be a two-body phase space multiplied by a polynomial function, which does not look to be easily written in one line, so I’ve chosen to do the following:

  • write a python function
  • make a TF1 out of it
  • turn it in a RooFit p.d.f.

I’m attaching a code example and a small dataset in the end of this post, but briefly it looks as the following:

import ROOT, math

...

def my_function(x, par):
    xx = x[0]
    a1 = par[0]
    a2 = par[1]
    ...
    f = f(x, a1, a2,...)
    return f

## name, function, left limit, right limit, number of parameters
f_ps = ROOT.TF1("f_ps", my_function, mass_left, mass_right, 3)

m_x = ROOT.RooRealVar("m_x", "m_x", mass_left, mass_right)

roo_ps_func = ROOT.RooFit.bindFunction(f_ps_pol2, m_x, pars)
roo_ps_pdf1 = ROOT.RooFit.bindPdf("roo_ps_pdf1",f_ps_pol2, m_x, a1, a2, a3)
roo_ps_pdf2 = ROOT.RooFit.bindPdf(f_ps_pol2, m_x)

where a1, a2 and a3 are RooRealVars for model parameters. I tried several methods but unfortunately none of them seems to allow fitting:

  • bound as a function
    roo_ps_func.fitTo(ds)
    leads to an error
    AttributeError: 'RooTFnBinding' object has no attribute 'fitTo'

  • bound as a pdf
    roo_ps_pdf1.fitTo(ds)
    leads to an error

File /snap/root-framework/954/usr/local/lib/ROOT/_pythonization/_roofit/_rooabspdf.py:64, in RooAbsPdf.fitTo(self, *args, **kwargs)
    60 The RooAbsPdf::fitTo() function is pythonized with the command argument pythonization.
    61 The keywords must correspond to the CmdArgs of the function.
    62 
    63 # Redefinition of `RooAbsPdf.fitTo` for keyword arguments.
--> 64 return self._fitTo["RooLinkedList const&"](args[0], _pack_cmd_args(*args[1:], **kwargs))

TypeError: Could not find "fitTo<RooLinkedList const&>" (set cppyy.set_debug() for C++ errors):
  Template method resolution failed:
  math domain error
  math domain error
  Failed to instantiate "operator()(double,double,double,double)"
  • other way of binding pdf
    roo_ps_pdf2.fitTo(ds)
    does not seem to see the parameters:
[#1] INFO:Fitting -- RooAbsPdf::fitTo(f_ps_pol2_over_f_ps_pol2_Int[m_x]) fixing normalization set for coefficient determination to observables in data
[#1] INFO:Fitting -- Creation of NLL object took 9.0357 ms
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_f_ps_pol2_over_f_ps_pol2_Int[m_x]_background) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#0] ERROR:Minimization -- RooMinimizer::fitFCN(): FCN function has zero parameters
[#0] ERROR:Minimization -- RooMinimizer: all function calls during minimization gave invalid NLL values!
[#0] WARNING:Minimization -- RooMinimizer::hesse: Error, run Migrad before Hesse!
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

So, it seems I’m missing something here. Please, could you give a hint of how to write a model to fit a distribution? Any other workaround would also be very appriciated.

If it is a necessary information, I’m using ROOT 6.36.04.

pdf_from_tf1.py (3.1 KB)
ds_bkg.root (59.8 KB)

Cheers,
Dasha.

2

Hello @meNikkie,

I will invite @jonas for this.

3

Dear @meNikkie,

before I answer to some of the Python specific problems: please try to avoid Python in the evaluation path of your model as much as you can. Especially if you wrap the Python function as a pdf that you minimize: there will be countless evaluations by the minimizer and for numeric integrals, so please try to stick C++ where you can, to get reasonable fit performance.

That being said, the best way to wrap custom C++ functions is RooGenericPdf for pdfs, and RooFormulaVar for functions. In your case, you should use the RooGenericPdf, and then fit it to the data.

Here is the full example:

import ROOT

m1 = m2 = 0.983 # protons

# range
mleft = m1 + m2
mright = m1 + m2 + 0.3

ROOT.gInterpreter.Declare(
    r"""
double phaseSpacePol2(double x, double a1, double a2, double a3)
{
    const double m1 = 0.983;
    const double m2 = 0.983;

    // phase space numerator: sqrt of the triangle function
    const double lambda = x*x*x*x + m1*m1*m1*m1 + m2*m2*m2*m2
                        - 2.*x*x*m1*m1 - 2.*m1*m1*m2*m2 - 2.*m2*m2*x*x;
    const double numr = TMath::Pi() * std::sqrt(lambda);

    // phase space denominator
    const double denr = 2. * x * x;

    // second order polynomial
    const double pol2 = a1*x*x + a2*x + a3;

    return pol2 * numr / denr;
}
"""
)

## turn into RooFit
# mass var
m_x  = ROOT.RooRealVar("m_x", "m_x", mleft, mright)

## 2nd order polynome coefficients with random values
a1 = ROOT.RooRealVar("a1" , "a1" , 3, -100, 100)
a2 = ROOT.RooRealVar("a2" , "a2" , 4, -100, 100)
a3 = ROOT.RooRealVar("a3" , "a3" , 5, -100, 100)
pars = ROOT.RooArgSet(a1, a2, a3)

# The PDF is just the C++ function referenced by name in the formula string.
# RooGenericPdf takes care of the normalization for you (it integrates the
# expression numerically over m_x), so you don't have to normalize by hand.
roo_ps_pdf = ROOT.RooGenericPdf(
    "roo_ps_pdf",
    "phaseSpacePol2(m_x, a1, a2, a3)",
    [m_x, a1, a2, a3],
)

## get dataset and try to fit
fil = ROOT.TFile('ds_bkg.root')
ds = fil.Get('background')

roo_ps_pdf.fitTo(ds)

So essentially you just declare the C++ function to the ROOT interpreter, and then you can readily use it in RooFit via RooGenericPdf.

But here is the actual problem in your example: the model is not valid. You give negative values to your sqrt() operation, which C++ will swallow and give you NaN (you fit won’t converge), and Python’s math module gives you this math domain error. If you don’t fix your model to be well defined for your orservable and parameter ranges, the fit won’t work. But to fix this is up to you, as you know best what measurement you’re doing :slightly_smiling_face: But I hope this clarifies how you define your own p.d.f. in RooFit.

Cheers,
Jonas

4

Dear @jonas

thank you very much for the explanation and especially for spending your time to create a working example!

Indeed, there was a typo in my script (proton mass of 983 MeV instead of 938), which shifted the fit range and lead to the errors. Correcting this everything works perfectly, thank you once again!

Cheers,
Dasha.