Skip to content

Commit 588ae6b

Browse files
committed
[RF] Don't double-scale pdfs plotted by RooSimultaneous::plotOn()
After a8ef8b0, when plotting an extended pdf it automatically scales itself to the number of expected events. However, when plotting slices of a RooSimultaneous, the normalization is already precomputed in `RooSimultanous::plotOn()` and should not be overridden again. To prevent this, flag the calculated normatization as the final number of events, and not some relative factor, which it is by default. A unit test that covers this case was also implemented. Closes #20383. (cherry picked from commit c948625)
1 parent 9b634c0 commit 588ae6b

File tree

3 files changed

+171
-5
lines changed

3 files changed

+171
-5
lines changed

roofit/roofitcore/src/RooSimultaneous.cxx

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -777,7 +777,8 @@ RooPlot* RooSimultaneous::plotOn(RooPlot *frame, RooLinkedList& cmdList) const
777777
projDataTmp(projData->reduce(RooFit::SelectVars(projDataVars), RooFit::Cut(cutString.c_str())));
778778

779779
// Override normalization and projection dataset
780-
RooCmdArg tmp1 = RooFit::Normalization(scaleFactor*wTable->getFrac(idxCatClone->getCurrentLabel()),stype) ;
780+
RooCmdArg tmp1 =
781+
RooFit::Normalization(scaleFactor * wTable->get(idxCatClone->getCurrentLabel()), RooAbsReal::NumEvent);
781782
RooCmdArg tmp2 = RooFit::ProjWData(*projDataSet,*projDataTmp) ;
782783

783784
// WVE -- do not adjust normalization for asymmetry plots

roofit/roofitcore/test/gtest_wrapper.h

Lines changed: 23 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,28 @@
3232

3333
#define ROOFIT_EVAL_BACKENDS ROOFIT_EVAL_BACKEND_LEGACY ROOFIT_EVAL_BACKEND_CUDA RooFit::EvalBackend::Cpu()
3434

35-
#define ROOFIT_EVAL_BACKENDS_WITH_CODEGEN ROOFIT_EVAL_BACKENDS, ROOFIT_EVAL_BACKEND_CODEGEN RooFit::EvalBackend::CodegenNoGrad()
35+
#define ROOFIT_EVAL_BACKENDS_WITH_CODEGEN \
36+
ROOFIT_EVAL_BACKENDS, ROOFIT_EVAL_BACKEND_CODEGEN RooFit::EvalBackend::CodegenNoGrad()
37+
38+
#include <gtest/gtest.h>
39+
#include <gmock/gmock.h>
40+
41+
#include <cmath>
42+
43+
MATCHER_P2(RelativeNear, expected, rel_tol,
44+
"is within relative tolerance around ref=" + ::testing::PrintToString(expected) +
45+
" (tol=" + ::testing::PrintToString(rel_tol) + ")")
46+
{
47+
const double diff = std::fabs(arg - expected);
48+
const double scale = std::fabs(expected);
49+
// One could also also consider the target value as scale, or an adaptive scale:
50+
// const double scale = std::max(std::fabs(arg), std::fabs(expected));
51+
52+
if (diff <= rel_tol * scale)
53+
return true;
54+
55+
*result_listener << "error relative to ref = " << diff / scale << ", absolute diff = " << diff;
56+
return false;
57+
}
3658

3759
#endif // RooFit_gtest_wrapper_h

roofit/roofitcore/test/testRooSimultaneous.cxx

Lines changed: 146 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -5,10 +5,12 @@
55
#include <RooAddPdf.h>
66
#include <RooAddition.h>
77
#include <RooCategory.h>
8+
#include <RooChebychev.h>
89
#include <RooConstVar.h>
910
#include <RooDataSet.h>
1011
#include <RooExponential.h>
1112
#include <RooFitResult.h>
13+
#include <RooGaussian.h>
1214
#include <RooGenericPdf.h>
1315
#include <RooHelpers.h>
1416
#include <RooMinimizer.h>
@@ -417,6 +419,8 @@ TEST(RooSimultaneous, ConditionalProdPdf)
417419
// channels can be extended. Also check if the likelihood can be created.
418420
TEST(RooSimultaneous, PartiallyExtendedPdfs)
419421
{
422+
RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
423+
420424
RooWorkspace ws;
421425
ws.factory("Gaussian::pdfA(x_a[-10, 10], mu_a[0, -10, 10], sigma_a[2.0, 0.1, 10.0])");
422426
ws.factory("Gaussian::pdfB(x_b[-10, 10], mu_b[0, -10, 10], sigma_b[2.0, 0.1, 10.0])");
@@ -427,7 +431,6 @@ TEST(RooSimultaneous, PartiallyExtendedPdfs)
427431
RooArgSet observables{*ws.var("x_a"), *ws.var("x_b"), *ws.cat("cat")};
428432

429433
auto &simPdf = *ws.pdf("simPdf");
430-
std::cout << simPdf.getVal() << std::endl;
431434

432435
// A completely extended pdf, just to easily create a toy dataset
433436
ws.factory("ExtendPdf::pdfAext(pdfA, n_b[1000., 100., 10000.])");
@@ -518,6 +521,7 @@ TEST_P(TestStatisticTest, RooSimultaneousSingleChannelCrossCheckWithCondVar)
518521
/// RooSimultaneous with the new CPU backend.
519522
TEST(RooSimultaneous, RangedExtendedRooAddPdf)
520523
{
524+
RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
521525

522526
const double nBkgA_nom = 9000;
523527
const double nBkgB_nom = 10000;
@@ -549,11 +553,12 @@ TEST(RooSimultaneous, RangedExtendedRooAddPdf)
549553
simPdf.getParameters(combData->get(), params);
550554
params.snapshot(paramsSnap);
551555

552-
Res fitResult{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(RooFit::EvalBackend::Cpu()))};
556+
Res fitResult{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Cpu()), PrintLevel(-1))};
553557

554558
params.assign(paramsSnap);
555559

556-
Res fitResultRef{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(RooFit::EvalBackend::Legacy()))};
560+
Res fitResultRef{
561+
simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Legacy()), PrintLevel(-1))};
557562

558563
EXPECT_TRUE(fitResult->isIdentical(*fitResultRef));
559564
}
@@ -563,6 +568,8 @@ TEST(RooSimultaneous, RangedExtendedRooAddPdf)
563568
/// with a projection dataset.
564569
TEST(RooSimultaneous, PlotProjWData)
565570
{
571+
RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
572+
566573
RooRealVar x("x", "x", -8, 8);
567574
x.setBins(1);
568575

@@ -590,11 +597,147 @@ TEST(RooSimultaneous, PlotProjWData)
590597
EXPECT_DOUBLE_EQ(frame->getCurve()->interpolate(0.), combData.sumEntries());
591598
}
592599

600+
/// Second part of GitHub issue #20383.
601+
/// Check that the the simultaneous pdf is normalized correctly to the data
602+
/// when plotting with a projection dataset, in the extended and non-extended
603+
/// case, based on the reproducer provided by the user who opened the issue.
604+
TEST(RooSimultaneous, PlotProjWDataExtended)
605+
{
606+
using namespace RooFit;
607+
608+
RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
609+
610+
RooRealVar xvar1("x1", "", 0.0, 10.0);
611+
RooRealVar xvar2("x2", "", 0.0, 10.0);
612+
613+
RooRealVar mean1("mean1", "", 3.0, 1.0, 4.0);
614+
RooRealVar mean2("mean2", "", 5.0, 4.0, 6.0);
615+
RooRealVar sigma("sigma", "", 1.0, 0.01, 2.0);
616+
617+
RooGaussian gauss1("gauss1", "", xvar1, mean1, sigma);
618+
RooGaussian gauss2("gauss2", "", xvar2, mean2, sigma);
619+
620+
RooRealVar a0("a0", "a0", -0.1, -1.0, 1.0);
621+
RooChebychev bkg1("bkg1", "b1", xvar1, {a0});
622+
623+
RooRealVar c0("c0", "c0", -0.1, -1.0, 1.0);
624+
RooChebychev bkg2("bkg2", "b2", xvar2, {c0});
625+
626+
RooRealVar s1("s1", "s1", 75.0, 0.0, 100000.0);
627+
RooRealVar b1("b1", "b1", 25.0, 0.0, 100000.0);
628+
629+
// Extended models
630+
RooAddPdf model1e("model1e", "", {gauss1, bkg1}, {s1, b1});
631+
632+
RooRealVar s2("s2", "s2", 50.0, 0.0, 100000.0);
633+
RooRealVar b2("b2", "b2", 50.0, 0.0, 100000.0);
634+
635+
RooAddPdf model2e("model2e", "model2e", {gauss2, bkg2}, {s2, b2});
636+
637+
// Non-extended models
638+
RooRealVar f1("f1", "f1", 0.75, 0.0, 1.0);
639+
RooAddPdf model1n("model1n", "model1n", {gauss1, bkg1}, {f1});
640+
641+
RooRealVar f2("f2", "f2", 0.50, 0.0, 1.0);
642+
RooAddPdf model2n("model2n", "model2n", RooArgList(gauss2, bkg2), RooArgList(f2));
643+
644+
// Case handling
645+
enum class Case {
646+
Gaussian,
647+
NonExtended,
648+
Extended
649+
};
650+
651+
auto caseName = [](Case c) {
652+
switch (c) {
653+
case Case::Gaussian: return "Gaussian";
654+
case Case::NonExtended: return "NonExtended";
655+
case Case::Extended: return "Extended";
656+
}
657+
return "Unknown";
658+
};
659+
660+
const std::vector<Case> cases = {Case::Gaussian, Case::NonExtended, Case::Extended};
661+
662+
// Helpers
663+
auto integrateLastCurve = [](RooPlot *plot) {
664+
const double xmin = plot->getPlotVar()->getMin();
665+
const double xmax = plot->getPlotVar()->getMax();
666+
return plot->getCurve()->average(xmin, xmax) * plot->getPlotVar()->numBins();
667+
};
668+
669+
constexpr double tol = 0.01; // tolerate 1 % sampling error
670+
671+
// Test body
672+
auto runCase = [&](Case c) {
673+
SCOPED_TRACE(std::string("Case = ") + caseName(c));
674+
675+
RooAbsPdf *model1 = nullptr;
676+
RooAbsPdf *model2 = nullptr;
677+
678+
switch (c) {
679+
case Case::Gaussian:
680+
model1 = &gauss1;
681+
model2 = &gauss2;
682+
break;
683+
case Case::NonExtended:
684+
model1 = &model1n;
685+
model2 = &model2n;
686+
break;
687+
case Case::Extended:
688+
model1 = &model1e;
689+
model2 = &model2e;
690+
break;
691+
}
692+
693+
ASSERT_NE(model1, nullptr);
694+
ASSERT_NE(model2, nullptr);
695+
696+
// Generate data
697+
std::unique_ptr<RooDataSet> data1{model1->generate(xvar1, 10000)};
698+
std::unique_ptr<RooDataSet> data2{model2->generate(xvar2, 1000)};
699+
700+
// Category
701+
RooCategory sample("sample", "");
702+
sample.defineType("Fit1");
703+
sample.defineType("Fit2");
704+
705+
RooDataSet data("combinedData", "", {xvar1, xvar2}, Index(sample),
706+
Import({{"Fit1", data1.get()}, {"Fit2", data2.get()}}));
707+
708+
// Simultaneous PDF
709+
RooSimultaneous sim_pdf("sim_pdf", "", sample);
710+
sim_pdf.addPdf(*model1, "Fit1");
711+
sim_pdf.addPdf(*model2, "Fit2");
712+
713+
std::unique_ptr<RooFitResult> result{sim_pdf.fitTo(data, Save(true), PrintLevel(-1))};
714+
715+
// Plot + checks
716+
RooPlot *frame1 = xvar1.frame();
717+
data.plotOn(frame1, Cut("sample==sample::Fit1"));
718+
sim_pdf.plotOn(frame1, Slice(sample, "Fit1"), ProjWData(sample, data));
719+
720+
EXPECT_THAT(integrateLastCurve(frame1), RelativeNear(data1->sumEntries(), tol));
721+
722+
RooPlot *frame2 = xvar2.frame();
723+
data.plotOn(frame2, Cut("sample==sample::Fit2"));
724+
sim_pdf.plotOn(frame2, Slice(sample, "Fit2"), ProjWData(sample, data));
725+
726+
EXPECT_THAT(integrateLastCurve(frame2), RelativeNear(data2->sumEntries(), tol));
727+
};
728+
729+
// Execute
730+
for (Case c : cases) {
731+
runCase(c);
732+
}
733+
}
734+
593735
/// JIRA ticket https://its.cern.ch/jira/browse/ROOT-7499
594736
/// Check that we can also generate Asimov datasets with non-integer weights
595737
/// via RooSimultaneous.
596738
TEST(RooSimultaneous, ExpectedDataWithNonIntegerWeights)
597739
{
740+
RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
598741

599742
RooWorkspace ws{"ws"};
600743
ws.factory("dummy_obs_a[0,1]");

0 commit comments

Comments
 (0)