Possible bug in 2D multi-range extended fit with EvalBackend("cpu")
Hi,
I am observing a discrepancy between the cpu and legacy RooFit backends when performing an extended fit in multiple disjoint ranges in 2D.
The issue is reproducible with a completely standalone script, without any external input files.
I tested it with:
-
ROOT 6.36.06
-
ROOT 6.40.02
and I obtain exactly the same behaviour in both versions.
The full-range fit gives identical results with the two backends, while the multi-range sideband fit gives very different fitted yields.
Here is the minimal reproducer:
import ROOT
def main():
ROOT.gROOT.SetBatch(True)
ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.ERROR)
ROOT.RooRandom.randomGenerator().SetSeed(12345)
ws = ROOT.RooWorkspace("ws")
for expression in (
"x[-5,5]",
"y[-5,5]",
"Gaussian::gx(x,mu[0],sigma[1])",
"Gaussian::gy(y,mu,sigma)",
"PROD::g(gx,gy)",
"Uniform::u({x,y})",
"SUM::model(ng[8000,0,30000]*g,nu[2000,0,30000]*u)",
):
ws.factory(expression)
obs = ROOT.RooArgSet(ws.var("x"), ws.var("y"))
bounds = [(-5, -2), (-2, 2), (2, 5)]
ranges = []
for i, x_bounds in enumerate(bounds):
for j, y_bounds in enumerate(bounds):
if i == 1 and j == 1:
continue
name = f"sb{i}{j}"
ranges.append(name)
ws.var("x").setRange(name, *x_bounds)
ws.var("y").setRange(name, *y_bounds)
data = ws.pdf("model").generate(obs, 15000)
n_sideband = data.sumEntries("abs(x)>2 || abs(y)>2")
print(
f"ROOT {ROOT.gROOT.GetVersion()}: "
f"full={data.numEntries()}, SB={n_sideband}",
flush=True,
)
results = {}
for restricted in (False, True):
for backend in ("cpu", "legacy"):
# Fresh copy prevents a preceding fit from changing range settings.
fit_ws = ws.Clone()
model = fit_ws.pdf("model")
options = (
[ROOT.RooFit.Range(",".join(ranges))]
if restricted else []
)
result = model.fitTo(
data,
ROOT.RooFit.Extended(True),
ROOT.RooFit.Save(),
ROOT.RooFit.PrintLevel(-1),
ROOT.RooFit.EvalBackend(backend),
*options,
)
yields = [
fit_ws.var(n).getVal()
for n in ("ng", "nu")
]
fit_obs = ROOT.RooArgSet(
fit_ws.var("x"),
fit_ws.var("y"),
)
fractions = [
sum(
fit_ws.pdf(pdf).createIntegral(
fit_obs,
ROOT.RooFit.NormSet(fit_obs),
ROOT.RooFit.Range(r),
).getVal()
for r in ranges
)
for pdf in ("g", "u")
]
expected_sb = sum(
n * fraction
for n, fraction in zip(yields, fractions)
)
label = "sidebands" if restricted else "full"
print(
f"{label:9s} {backend:6s}: "
f"N_g={yields[0]:.3f}, "
f"N_u={yields[1]:.3f}, "
f"expected_SB={expected_sb:.3f}, "
f"status={result.status()}, "
f"covQual={result.covQual()}",
flush=True,
)
results[restricted, backend] = yields
if __name__ == "__main__":
main()
With ROOT 6.40.02 I get:
ROOT 6.40.02: full=15000, SB=3431.0
full cpu : N_g=12146.553, N_u=2854.253,
expected_SB=3477.756,
status=0, covQual=3
full legacy: N_g=12146.553, N_u=2854.253,
expected_SB=3477.756,
status=0, covQual=3
sidebands cpu : N_g=2818.327, N_u=729.358,
expected_SB=863.292,
status=0, covQual=3
sidebands legacy: N_g=11199.522, N_u=2899.005,
expected_SB=3431.129,
status=0, covQual=3
The full-range fits agree perfectly.
For the sideband fit, however, the legacy backend gives
expected_SB = 3431.129
which agrees with the actual number of events in the sidebands,
N_SB = 3431
while the cpu backend gives
expected_SB = 863.292.
Interestingly, the fitted yields differ essentially by a common scale factor:
11199.522 / 2818.327 ≈ 3.97
2899.005 / 729.358 ≈ 3.97
so the discrepancy looks more like a normalization / extended-term issue than a different likelihood minimum.
Both fits report
status = 0
covQual = 3
and the discrepancy is still present in ROOT 6.40.02.
Is this expected behaviour for Extended(True) combined with a comma-separated multi-range fit in more than one dimension, or is this a bug in the cpu backend normalization?
In my understanding, the fitted extended coefficients should correspond consistently to the full-range yields, with the expected number of events in the fitted sideband region obtained from the component integrals over the selected ranges.
Thanks!