Failing numeric integration of custom PDF under different ROOT versions

I have problems with a custom written PDF and its integration. These problems showed up when I switched from ROOT 5.34.18 to version 6.04.02 but they are also already present with version 5.34.30.
Nothing in the changelogs gave me a clue to why this problem appeared. I am using pyROOT.

My custom PDF is a 2d Gaussian with additional features (factors to scale the center in x and y, efficiencies in both variables, input number of events to make it extended). If I plot the PDF for certain positions of the peak, the numeric integration seems to fail and thus the normalization too. I can see this in the plotted PDF and the integral I get from the expectedEvents method. With ROOT version 5.34.18 there is no problem.

It would be great if someone has an idea of why this problem shows up and how I can fix it. I would really like to do the switch to ROOT 6 but for now I’m stuck with the old version.

Here are some simple plots from the script and example results to illustrate the problem:

ROOT 5.34.30
[color=#00BF00]OK when peak @ (x=2,y=2)
inputevents = 123
integral = 123.000021457[/color]

[color=#FF0000]Problems with peak = (x=10,y=10)
inputevents = 123
integral = 2.47409635955e-30[/color]

ROOT 5.34.18
[color=#00BF00]OK when peak = (x=2,y=2)
inputevents = 123
integral = 123.000021457[/color]

[color=#00BF00]OK when peak = (x=10,y=10)
inputevents = 123
integral = 122.99999739[/color]

I attached a minimal working example in a zip-file. It contains the short python script to plot the PDF and the C++ class with the PDF description:

mwe.zip (1.98 KB)

The script does nothing more than to create fake data and plot a 2D histogram as well as the 2 projections in x and y with the fake data.

[code]#!/usr/bin/env python

import ROOT
from ROOT import RooRealVar, RooArgSet, TCanvas
RF = ROOT.RooFit

peak = 12
inputevents = 123

x = RooRealVar(“x”, “x”, 0, 15)
y = RooRealVar(“y”, “y”, 0, 15)
observables = RooArgSet(x,y)

pdf = ROOT.TestPdf(“pdf”, “pdf”, x, y, RF.RooConst(peak), RF.RooConst(0.2), RF.RooConst(0.2), RF.RooConst(1.0), RF.RooConst(1.0), RF.RooConst(1.0), RF.RooConst(1.0), RF.RooConst(inputevents))

hist = pdf.createHistogram(“hist”, x, RF.YVar(y), RF.Extended(1))
data = pdf.generate(observables, 1000)
datahist = data.createHistogram(x, y, 1500, 1500)

c1 = TCanvas(“c1”, “c1”, 1200, 400)
c1.Divide(3,1)

c1.cd(1)
c1.cd(1).SetLogz()
hist.Draw(“COLZ”)
datahist.Draw(“PSAME”)

c1.cd(2)
xframe = x.frame()
data.plotOn(xframe)
pdf.plotOn(xframe)
xframe.Draw()

c1.cd(3)
yframe = y.frame()
data.plotOn(yframe)
pdf.plotOn(yframe)
yframe.Draw()

print “inputevents =”, inputevents
print “integral =”, pdf.expectedEvents(observables)
[/code]

Hi,

The problem is caused by the difficulty of the numerical integration method to integrate a so small peak in a very large range. If you can, you should reduce the integration range of your function.
You might try also to use the MC integrator method (I have tried and it seems to work in your case).
You can do this by adding the line:

``ROOT.RooAbsReal.defaultIntegratorConfig().method2D().setLabel("RooMCIntegrator")``

Best Regards

Lorenzo

Great, thank you! With that integration method it works, even with ROOT 6:

ROOT 6.04.02:
[color=#00BF40]OK when peak @ (x=10,y=10)
inputevents = 123
integral = 123.492598438[/color]

But the precision of this integration is much worse now than before. I am afraid this will not be sufficient for my purposes. The PDF I have here is part of a large RooAddPdf with ~ 20 components in 2 observables. Is there any other way I can keep sufficient precision without failing to integrate over a range like here?

Hi,

Yes, this is a problem with the MC integration. You can try to increase the number of function calls by changing this configuration option for the ROoMCIntegration. You need however to not be patient, the CPU time will increase as well.
For example for changing from the default value of 5000 to 50000 you do:

``ROOT.RooAbsReal.defaultIntegratorConfig().getConfigSection("RooMCIntegrator").setRealValue("nIntPerDim",50000)``

Lorenzo

Thanks, I played around with this setting but the results are not entirely convincing.

With 50,000 I get a 0.1/123 error:
inputevents = 123
integral = 122.900212711

With 500,000 it looks better and is still quite fast:
inputevents = 123
integral = 122.987608925

to reach the same precision as I had with the older ROOT version I need >1 million.

What I don’t understand is why everything worked well with the older ROOT version (without the MC integraton method I suppose). What can have changed?
Amd I am not sure how to proceed: Switch to ROOT 6 and use the MC integration with many samples (?) (I would have to set this manually in all my scripts, and make sure to never forget) or stay with ROOT 5.34.18 with which everything worked fine directly.

I could run 5.34.18 and in did it seems to work with that version. I will investigate it

Lorenzo

HI,

I found the cause of the problem. It is in having introduced the use of absolute tolerance in the AdaptiveIntegrator after 5.34.18
I will fix this for the next version, but you can use this workaround:

``ROOT.RooAbsReal.defaultIntegratorConfig().setEpsAbs(1E-300)``

You might want to change later to the default value (1.E-9) if you are using later other integrators
Thank you for having reported this important problem.

Best Regards

Lorenzo

Hi Lorenzo,

that’s great news! Thank you for the quick solution, I’ll use the workaround and switch to ROOT 6 without any worries then
Best,

Lukas