-
Notifications
You must be signed in to change notification settings - Fork 1.4k
Open
Labels
Description
Description
This was reported on the forum:
The reproducer below is a simplified version of the reproducer on the forum, translated to C++ so it can later be used as a unit test.
Reproducer
#include "RooRealVar.h"
#include "RooCategory.h"
#include "RooDataHist.h"
#include "RooArgSet.h"
#include "RooPolyVar.h"
#include "RooHistPdf.h"
#include "RooGaussian.h"
#include "RooFFTConvPdf.h"
#include "RooSimultaneous.h"
#include "TString.h"
#include <iostream>
#include <vector>
void ana_globalFit()
{
using namespace RooFit;
// create all RooRealVar
RooRealVar xvar("x", "", 0.0, 10.0);
// Define category to distinguish physics and control samples events
RooCategory sample("sample", "sample");
sample.defineType("sample_0");
// create all RooDataHist for DATA, prepare dictionary
RooDataHist data("data_0", "", RooArgList(xvar));
std::vector<int> content = {0, 0, 0, 0, 1, 0, 2, 0, 0, 3, 5, 8, 6, 15, 17, 18, 23, 14, 24, 25,
33, 33, 43, 39, 41, 48, 51, 39, 43, 39, 49, 27, 31, 29, 27, 29, 38, 34, 22, 23,
21, 12, 19, 12, 10, 10, 11, 9, 10, 9, 8, 6, 4, 7, 4, 6, 6, 3, 2, 2,
4, 4, 4, 5, 2, 2, 1, 1, 1, 1, 1, 0, 3, 1, 0, 2, 0, 1, 0, 0,
1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1};
for (size_t i = 0; i < content.size(); ++i) {
data.set(i, content[i], -1);
}
xvar.setBins(data.numEntries());
// create the GLOBAL dataset
std::map<std::string, std::unique_ptr<RooDataHist>> importMap;
importMap["range_0"] = std::make_unique<RooDataHist>(data);
RooDataHist dataALL("data", "", RooArgList(xvar), Index(sample), Import(importMap));
// create the RooRealVar for MC
RooRealVar E("E_0", "", 0, 500);
E.setBins(10000, "cache");
std::vector<int> content2 = {
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 13, 19, 55, 114,
145, 203, 295, 361, 497, 607, 683, 885, 964, 1167, 1240, 1328, 1379, 1479, 1583, 1638, 1667, 1667, 1705,
1694, 1677, 1757, 1723, 1556, 1606, 1563, 1513, 1516, 1346, 1329, 1336, 1322, 1252, 1170, 1124, 1092, 1060, 1017,
999, 959, 881, 836, 853, 772, 784, 747, 716, 654, 602, 602, 588, 535, 549, 517, 479, 495, 425,
402, 449, 403, 355, 348, 393, 361, 332, 320, 325, 277, 312, 277, 282, 230, 256, 253, 234, 234,
204, 205, 204, 195, 215, 194, 208, 183, 167, 144, 164, 159, 168, 153, 148, 161, 140, 135, 126,
149, 137, 115, 121, 110, 117, 104, 108, 121, 122, 83, 97, 91, 85, 81, 66, 66, 71, 80,
76, 75, 74, 73, 76, 53, 69, 64, 67, 70, 77, 76, 48, 55, 45, 56, 62, 50, 51,
49, 44, 53, 47, 48, 55, 49, 56, 41, 43, 33, 57, 53, 40, 44, 28, 42, 34, 35,
44, 45, 32, 34, 35, 28, 39, 27, 34, 21, 31, 26, 24, 28, 41, 31, 28, 30, 29,
30, 15, 23, 25, 27, 27, 21, 32, 20, 22, 32, 23, 23, 17, 20, 29, 23, 21, 21,
22, 18, 20, 27, 17, 20, 16, 18, 31, 20, 11, 15, 17, 20, 16, 11, 15, 14, 12,
17, 18, 13, 14, 13, 14, 14, 12, 13, 13, 9, 13, 12, 9, 7, 12, 17, 14, 6,
7, 17, 9, 10, 18, 11, 21, 10, 15, 8, 16, 11, 19, 16, 10, 13, 13, 11, 7,
14, 12, 12, 18, 11, 8, 9, 4, 9, 9, 4, 12, 8, 11, 7, 5, 7, 7, 7,
9, 7, 8, 6, 5, 5, 6, 11, 11, 11, 9, 10, 7, 6, 7, 8, 3, 8, 2,
7, 7, 10, 6, 8, 10, 9, 10, 11, 9, 7, 6, 15, 6, 7, 6, 5, 4, 8,
6, 7, 3, 2, 6, 4, 3, 10, 7, 4, 5, 6, 6, 5, 5, 2, 5, 6, 10,
1};
// create all RootDataHist for MC
RooRealVar Ehist("E_0", "", 0, 10);
Ehist.setBins(content2.size());
RooDataHist dh("dh_0", "", RooArgList(Ehist));
for (size_t i = 0; i < content2.size(); ++i) {
dh.set(i, content2[i], -1);
}
// create the scale variable for MC - assume a single scale is ok for all energy ranges.
RooRealVar scale("scale_0", "", 0.95, 0.5, 1.5);
RooRealVar p0("p0", "", 0.0);
RooPolyVar Qf("HCAL_F_0", "", xvar, RooArgSet(p0, scale));
// create all histo PDFs for MC.
RooHistPdf histpdf("histpdf_0", "", E, dh);
// gaussian model
RooRealVar mg("mg", "", 0);
RooRealVar sg("sg_0", "", 0.7, 0.1, 10.);
RooGaussian gauss("gauss_0", "", E, mg, sg);
// create all histo PDFs for MC after smearing.
RooFFTConvPdf lxg("lxg_0", "", Qf, E, histpdf, gauss);
// prepare the simultaneous PDF
RooSimultaneous simPDF("simPDF", "", sample);
simPDF.addPdf(lxg, "sample_0");
// compute NLLs
std::unique_ptr<RooAbsReal> nll_single{lxg.createNLL(data, EvalBackend("legacy"))};
std::unique_ptr<RooAbsReal> nll{simPDF.createNLL(dataALL, EvalBackend("legacy"))};
std::vector<double> test_vals = {0.2, 0.3};
for (auto val : test_vals) {
sg.setVal(val);
double val_single = nll_single->getVal();
double val_sim = nll->getVal();
std::cout << val_single << std::endl;
std::cout << val_sim << std::endl;
}
}The output shows that the simultaneous NLL is always zero, while it's expected to be the same:
------------------------------------------------------------------
| Welcome to ROOT 6.37.01 https://root.cern |
| (c) 1995-2025, The ROOT Team; conception: R. Brun, F. Rademakers |
| Built for linuxx8664gcc on Jan 01 1980, 00:00:00 |
| From heads/master@v6-37-01-8453-g301e0782d63 |
| With clang version 21.1.1 |
| Try '.help'/'.?', '.demo', '.license', '.credits', '.quit'/'.q' |
------------------------------------------------------------------
root [0]
Processing ana_globalFit.C...
[#1] INFO:InputArguments -- RooDataHist::importDHistSet(data) defining state "range_0" in index category sample
[#1] INFO:Eval -- RooRealVar::setRange(E_0) new range named 'refrange_fft_lxg_0' created with bounds [0,500]
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(lxg_0) creating new cache 0x5595b0492960 with pdf histpdf_0_CONV_gauss_0_CACHE_Obs[p0,scale_0,x]_NORM_x for nset (x) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(histpdf_0_CONV_gauss_0_CACHE_Obs[p0,scale_0,x]_NORM_x_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
[#1] INFO:Fitting -- Creation of NLL object took 39.9617 ms
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 0 slave calculators.
[#1] INFO:Fitting -- Creation of NLL object took 105.859 μs
2236.42
0
1927.05
0
root [1]The new default "cpu" backend has the same problem.
ROOT version
Probably any.