From 8087814dfc69db5e8950cd834c0bc88ae811d96f Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sat, 5 Sep 2026 19:49:28 +0000 Subject: [PATCH] [RF] Convert stressRooFit tests to reference-free gtests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Convert 15 stressRooFit tests to google tests that validate against analytic expectations, deterministic data, or self-consistency instead of the hardcoded stressRooFit_ref.root file, following the precedent of the rf705 conversion to testRooIntegralMorph. This way, expected numerical changes can be distinguished from real regressions. - rf110/rf308 (1D/2D normalization, createCdf), rf111 (numeric integration config) and rf314 (parameterized-range fit): testRooAbsPdf - rf302 (utility function composition): testRooAbsReal - rf402 (dataset reduce/merge/append): testRooDataSet - rf404-rf406 (Roo1DTable, category ranges, threshold/binning/mapped/ super/multi categories): testRooCategory - rf109 (RooPlot::chiSquare, residHist, pullHist): testRooHist - rf605 (RooProfileLL): new testRooProfileLL - rf701-rf703 (RooEfficiency 1D/2D, RooEffProd): new testRooEfficiency - rf704 (amplitude sum pdf): new testRooRealSumPdf Several of these classes had no google test coverage at all before (RooProfileLL, RooEfficiency, RooEffProd, the category mapping classes, RooPlot residual/pull machinery, and dataset merge/append). The fit-based conversions are parametrized over the RooFit evaluation backends via a new shared fixture in gtest_wrapper.h, matching the backend coverage of the removed stressRooFit invocations (including codegen for the RooEfficiency fits, which skips only the RooEffProd case that the stress suite also excluded). Tolerances for comparing results of independent minimizations are kept at the scale of the Minuit convergence criterion to avoid flakiness. The converted tests are removed from the stressRooFit suite; the corresponding entries in stressRooFit_ref.root are simply unused, so the reference file does not need to be regenerated. 🤖 Done with the help of AI --- roofit/roofitcore/test/CMakeLists.txt | 3 + roofit/roofitcore/test/gtest_wrapper.h | 40 + roofit/roofitcore/test/stressRooFit.cxx | 15 - roofit/roofitcore/test/stressRooFit_tests.h | 1174 ------------------ roofit/roofitcore/test/testRooAbsPdf.cxx | 170 +++ roofit/roofitcore/test/testRooAbsReal.cxx | 66 + roofit/roofitcore/test/testRooCategory.cxx | 254 ++++ roofit/roofitcore/test/testRooDataSet.cxx | 71 ++ roofit/roofitcore/test/testRooEfficiency.cxx | 160 +++ roofit/roofitcore/test/testRooHist.cxx | 72 ++ roofit/roofitcore/test/testRooProfileLL.cxx | 107 ++ roofit/roofitcore/test/testRooRealSumPdf.cxx | 94 ++ 12 files changed, 1037 insertions(+), 1189 deletions(-) create mode 100644 roofit/roofitcore/test/testRooEfficiency.cxx create mode 100644 roofit/roofitcore/test/testRooProfileLL.cxx create mode 100644 roofit/roofitcore/test/testRooRealSumPdf.cxx diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index ac9d5bb73b7f0..810eb310a7fb4 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -35,6 +35,7 @@ ROOT_ADD_GTEST(testRooFitCore testRooLinkedList.cxx testRooPolyFunc.cxx testRooProdPdf.cxx + testRooRealSumPdf.cxx testRooSTLRefCountList.cxx testRooWrapperPdf.cxx testSimple.cxx @@ -81,10 +82,12 @@ if(clad) ROOT_ADD_GTEST(testRooFuncWrapper testRooFuncWrapper.cxx LIBRARIES RooFitCore RooFit HistFactory) endif() ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testRooEfficiency testRooEfficiency.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooExtendedBinding testRooExtendedBinding.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooMCStudy testRooMCStudy.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooMinimizer testRooMinimizer.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooMulti testRooMulti.cxx LIBRARIES RooFitCore RooFit) +ROOT_ADD_GTEST(testRooProfileLL testRooProfileLL.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooRombergIntegrator testRooRombergIntegrator.cxx LIBRARIES MathCore RooFitCore) ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooTruthModel testRooTruthModel.cxx LIBRARIES RooFitCore RooFit diff --git a/roofit/roofitcore/test/gtest_wrapper.h b/roofit/roofitcore/test/gtest_wrapper.h index f342c55575ab8..7a62f5bc93530 100644 --- a/roofit/roofitcore/test/gtest_wrapper.h +++ b/roofit/roofitcore/test/gtest_wrapper.h @@ -35,10 +35,50 @@ #define ROOFIT_EVAL_BACKENDS_WITH_CODEGEN \ ROOFIT_EVAL_BACKENDS, ROOFIT_EVAL_BACKEND_CODEGEN RooFit::EvalBackend::CodegenNoGrad() +#include +#include +#include +#include + #include #include #include +#include +#include + +/// Common fixture for tests that are parametrized over the RooFit evaluation +/// backends: fixes the random seed and silences RooFit messages below WARNING. +class RooFitEvalBackendTest : public testing::TestWithParam> { +public: + RooFitEvalBackendTest() : _evalBackend{RooFit::EvalBackend::Legacy()} {} + +private: + void SetUp() override + { + RooRandom::randomGenerator()->SetSeed(1337ul); + _evalBackend = std::get<0>(GetParam()); + _changeMsgLvl = std::make_unique(RooFit::WARNING); + } + + void TearDown() override { _changeMsgLvl.reset(); } + +protected: + RooFit::EvalBackend _evalBackend; + +private: + std::unique_ptr _changeMsgLvl; +}; + +/// Check that the floating fit parameter with the given name is within +/// nSigma fit errors of the truth value. +inline void expectParamNear(RooFitResult const &res, const char *name, double truthVal, double nSigma = 5.) +{ + auto *param = static_cast(res.floatParsFinal().find(name)); + ASSERT_NE(param, nullptr) << name; + EXPECT_GT(param->getError(), 0.) << name; + EXPECT_NEAR(param->getVal(), truthVal, nSigma * param->getError()) << name; +} MATCHER_P2(RelativeNear, expected, rel_tol, "is within relative tolerance around ref=" + ::testing::PrintToString(expected) + diff --git a/roofit/roofitcore/test/stressRooFit.cxx b/roofit/roofitcore/test/stressRooFit.cxx index 17bcfb842b589..510cd95534172 100644 --- a/roofit/roofitcore/test/stressRooFit.cxx +++ b/roofit/roofitcore/test/stressRooFit.cxx @@ -106,9 +106,6 @@ int stressRooFit(const char *refFile, bool writeRef, int doVerbose, int oneTest, testList.push_back(new TestBasic103(fref, writeRef, doVerbose)); testList.push_back(new TestBasic105(fref, writeRef, doVerbose)); testList.push_back(new TestBasic108(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic109(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic110(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic111(fref, writeRef, doVerbose)); testList.push_back(new TestBasic201(fref, writeRef, doVerbose)); testList.push_back(new TestBasic202(fref, writeRef, doVerbose)); testList.push_back(new TestBasic203(fref, writeRef, doVerbose)); @@ -117,36 +114,24 @@ int stressRooFit(const char *refFile, bool writeRef, int doVerbose, int oneTest, testList.push_back(new TestBasic208(fref, writeRef, doVerbose)); testList.push_back(new TestBasic209(fref, writeRef, doVerbose)); testList.push_back(new TestBasic301(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic302(fref, writeRef, doVerbose)); testList.push_back(new TestBasic303(fref, writeRef, doVerbose)); testList.push_back(new TestBasic304(fref, writeRef, doVerbose)); testList.push_back(new TestBasic305(fref, writeRef, doVerbose)); testList.push_back(new TestBasic306(fref, writeRef, doVerbose)); testList.push_back(new TestBasic307(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic308(fref, writeRef, doVerbose)); testList.push_back(new TestBasic310(fref, writeRef, doVerbose)); testList.push_back(new TestBasic311(fref, writeRef, doVerbose)); testList.push_back(new TestBasic312(fref, writeRef, doVerbose)); testList.push_back(new TestBasic313(fref, writeRef, doVerbose)); testList.push_back(new TestBasic315(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic314(fref, writeRef, doVerbose)); testList.push_back(new TestBasic316(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic402(fref, writeRef, doVerbose)); testList.push_back(new TestBasic403(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic404(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic405(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic406(fref, writeRef, doVerbose)); testList.push_back(new TestBasic501(fref, writeRef, doVerbose)); testList.push_back(new TestBasic599(fref, writeRef, doVerbose)); testList.push_back(new TestBasic602(fref, writeRef, doVerbose)); testList.push_back(new TestBasic604(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic605(fref, writeRef, doVerbose)); testList.push_back(new TestBasic606(fref, writeRef, doVerbose)); testList.push_back(new TestBasic607(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic701(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic702(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic703(fref, writeRef, doVerbose)); - testList.push_back(new TestBasic704(fref, writeRef, doVerbose)); testList.push_back(new TestBasic706(fref, writeRef, doVerbose)); testList.push_back(new TestBasic707(fref, writeRef, doVerbose)); testList.push_back(new TestBasic708(fref, writeRef, doVerbose)); diff --git a/roofit/roofitcore/test/stressRooFit_tests.h b/roofit/roofitcore/test/stressRooFit_tests.h index 71d4b5258fb66..d43acce667d9d 100644 --- a/roofit/roofitcore/test/stressRooFit_tests.h +++ b/roofit/roofitcore/test/stressRooFit_tests.h @@ -1,4 +1,3 @@ -#include #include #include #include @@ -9,7 +8,6 @@ #include #include #include -#include #include #include #include @@ -20,8 +18,6 @@ #include #include #include -#include -#include #include #include #include @@ -37,9 +33,7 @@ #include #include #include -#include #include -#include #include #include #include @@ -47,15 +41,10 @@ #include #include #include -#include #include -#include -#include #include #include -#include #include -#include #include #include #include @@ -539,216 +528,6 @@ class TestBasic108 : public RooUnitTest { } }; -// Calculating chi^2 from histograms and curves in RooPlots, making histogram of residual and pull distributions -class TestBasic109 : public RooUnitTest { -public: - TestBasic109(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Calculation of chi^2 and residuals in plots", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // S e t u p m o d e l - // --------------------- - - // Create observables - RooRealVar x("x", "x", -10, 10); - - // Create Gaussian - RooRealVar sigma("sigma", "sigma", 3, 0.1, 10); - RooRealVar mean("mean", "mean", 0, -10, 10); - RooGaussian gauss("gauss", "gauss", x, 0.0, sigma); - - // Generate a sample of 1000 events with sigma=3 - std::unique_ptr data{gauss.generate(x, 10000)}; - - // Change sigma to 3.15 - sigma = 3.15; - - // P l o t d a t a a n d s l i g h t l y d i s t o r t e d m o d e l - // --------------------------------------------------------------------------- - - // Overlay projection of gauss with sigma=3.15 on data with sigma=3.0 - RooPlot *frame1 = x.frame(Title("Data with distorted Gaussian pdf"), Bins(40)); - data->plotOn(frame1, DataError(RooAbsData::SumW2)); - gauss.plotOn(frame1); - - // C a l c u l a t e c h i ^ 2 - // ------------------------------ - - // Show the chi^2 of the curve w.r.t. the histogram - // If multiple curves or datasets live in the frame you can specify - // the name of the relevant curve and/or dataset in chiSquare() - regValue(frame1->chiSquare(), "rf109_chi2"); - - // S h o w r e s i d u a l a n d p u l l d i s t s - // ------------------------------------------------------- - - // Construct a histogram with the residuals of the data w.r.t. the curve - // we set `useAverage` to false for this test because this was done for the reference histogram - RooHist *hresid = frame1->residHist(nullptr, nullptr, false, false); - - // Construct a histogram with the pulls of the data w.r.t the curve - // we set `useAverage` to false for this test because this was done for the reference histogram - RooHist *hpull = frame1->pullHist(nullptr, nullptr, false); - - // Create a new frame to draw the residual distribution and add the distribution to the frame - RooPlot *frame2 = x.frame(Title("Residual Distribution")); - frame2->addPlotable(hresid, "P"); - - // Create a new frame to draw the pull distribution and add the distribution to the frame - RooPlot *frame3 = x.frame(Title("Pull Distribution")); - frame3->addPlotable(hpull, "P"); - - regPlot(frame1, "rf109_plot1"); - regPlot(frame2, "rf109_plot2"); - regPlot(frame3, "rf109_plot3"); - - // delete hresid ; - // delete hpull ; - - return true; - } -}; - -// Examples on normalization of p.d.f.s, integration of p.d.fs, construction of -// cumulative distribution functions from p.d.f.s in one dimension. -class TestBasic110 : public RooUnitTest { -public: - TestBasic110(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Normalization of p.d.f.s in 1D", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // S e t u p m o d e l - // --------------------- - - // Create observables x,y - RooRealVar x("x", "x", -10, 10); - - // Create p.d.f. gaussx(x,-2,3) - RooGaussian gx("gx", "gx", x, -2.0, 3.0); - - // R e t r i e v e r a w & n o r m a l i z e d v a l u e s o f R o o F i t p . d . f . s - // -------------------------------------------------------------------------------------------------- - - // Return 'raw' unnormalized value of gx - regValue(gx.getVal(), "rf110_gx"); - - // Return value of gx normalized over x in range [-10,10] - RooArgSet nset(x); - - regValue(gx.getVal(&nset), "rf110_gx_Norm[x]"); - - // Create object representing integral over gx - // which is used to calculate gx_Norm[x] == gx / gx_Int[x] - std::unique_ptr igx{gx.createIntegral(x)}; - regValue(igx->getVal(), "rf110_gx_Int[x]"); - - // I n t e g r a t e n o r m a l i z e d p d f o v e r s u b r a n g e - // ---------------------------------------------------------------------------- - - // Define a range named "signal" in x from -5,5 - x.setRange("signal", -5, 5); - - // Create an integral of gx_Norm[x] over x in range "signal" - // This is the fraction of of p.d.f. gx_Norm[x] which is in the - // range named "signal" - std::unique_ptr igx_sig{gx.createIntegral(x, NormSet(x), Range("signal"))}; - regValue(igx_sig->getVal(), "rf110_gx_Int[x|signal]_Norm[x]"); - - // C o n s t r u c t c u m u l a t i v e d i s t r i b u t i o n f u n c t i o n f r o m p d f - // ----------------------------------------------------------------------------------------------------- - - // Create the cumulative distribution function of gx - // i.e. calculate Int[-10,x] gx(x') dx' - std::unique_ptr gx_cdf{gx.createCdf(x)}; - - // Plot cdf of gx versus x - RooPlot *frame = x.frame(Title("c.d.f of Gaussian p.d.f")); - gx_cdf->plotOn(frame); - - regPlot(frame, "rf110_plot1"); - - return true; - } -}; - -// Configuration and customization of how numeric (partial) integrals -// are executed. -class TestBasic111 : public RooUnitTest { -public: - TestBasic111(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Numeric integration configuration", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // A d j u s t g l o b a l 1 D i n t e g r a t i o n p r e c i s i o n - // ---------------------------------------------------------------------------- - - // Example: Change global precision for 1D integrals from 1e-7 to 1e-6 - // - // The relative epsilon (change as fraction of current best integral estimate) and - // absolute epsilon (absolute change w.r.t last best integral estimate) can be specified - // separately. For most p.d.f integrals the relative change criterium is the most important, - // however for certain non-p.d.f functions that integrate out to zero a separate absolute - // change criterium is necessary to declare convergence of the integral - // - // NB: This change is for illustration only. In general the precision should be at least 1e-7 - // for normalization integrals for MINUIT to succeed. - // - RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1e-6); - RooAbsReal::defaultIntegratorConfig()->setEpsRel(1e-6); - - // N u m e r i c i n t e g r a t i o n o f l a n d a u p d f - // ------------------------------------------------------------------ - - RooRealVar x("x", "x", -10, 10); - RooLandau landau("landau", "landau", x, 0.0, 0.1); - - // Disable analytic integration from demonstration purposes - landau.forceNumInt(true); - - // Calculate integral over landau with default choice of numeric integrator - std::unique_ptr intLandau{landau.createIntegral(x)}; - double val = intLandau->getVal(); - regValue(val, "rf111_val1"); - - // S a m e w i t h c u s t o m c o n f i g u r a t i o n - // ----------------------------------------------------------- - - // Construct a custom configuration which uses the adaptive Gauss-Kronrod technique - // for closed 1D integrals - RooNumIntConfig customConfig(*RooAbsReal::defaultIntegratorConfig()); -#ifdef ROOFITMORE - customConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - - // Calculate integral over landau with custom integral specification - std::unique_ptr intLandau2{landau.createIntegral(x, NumIntConfig(customConfig))}; - double val2 = intLandau2->getVal(); - regValue(val2, "rf111_val2"); - - // A d j u s t i n g d e f a u l t c o n f i g f o r a s p e c i f i c p d f - // ------------------------------------------------------------------------------------- - - // Another possibility: associate custom numeric integration configuration as default for object 'landau' - landau.setIntegratorConfig(customConfig); - - // Calculate integral over landau custom numeric integrator specified as object default - std::unique_ptr intLandau3{landau.createIntegral(x)}; - double val3 = intLandau3->getVal(); - regValue(val3, "rf111_val3"); - - return true; - } -}; - // Composite p.d.f with signal and background component. class TestBasic201 : public RooUnitTest { public: @@ -1372,87 +1151,6 @@ class TestBasic301 : public RooUnitTest { } }; -// Utility functions classes available for use in tailoring -// of composite (multidimensional) pdfs -class TestBasic302 : public RooUnitTest { -public: - TestBasic302(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Sum and product utility functions", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // C r e a t e o b s e r v a b l e s , p a r a m e t e r s - // ----------------------------------------------------------- - - // Create observables - RooRealVar x("x", "x", -5, 5); - RooRealVar y("y", "y", -5, 5); - - // Create parameters - RooRealVar a0("a0", "a0", -1.5, -5, 5); - RooRealVar a1("a1", "a1", -0.5, -1, 1); - RooRealVar sigma("sigma", "width of gaussian", 0.5); - - // U s i n g R o o F o r m u l a V a r t o t a i l o r p d f - // ----------------------------------------------------------------------- - - // Create interpreted function f(y) = a0 - a1*sqrt(10*abs(y)) - RooFormulaVar fy_1("fy_1", "a0-a1*sqrt(10*abs(y))", RooArgSet(y, a0, a1)); - - // Create gauss(x,f(y),s) - RooGaussian model_1("model_1", "Gaussian with shifting mean", x, fy_1, sigma); - - // U s i n g R o o P o l y V a r t o t a i l o r p d f - // ----------------------------------------------------------------------- - - // Create polynomial function f(y) = a0 + a1*y - RooPolyVar fy_2("fy_2", "fy_2", y, RooArgSet(a0, a1)); - - // Create gauss(x,f(y),s) - RooGaussian model_2("model_2", "Gaussian with shifting mean", x, fy_2, sigma); - - // U s i n g R o o A d d i t i o n t o t a i l o r p d f - // ----------------------------------------------------------------------- - - // Create sum function f(y) = a0 + y - RooAddition fy_3("fy_3", "a0+y", RooArgSet(a0, y)); - - // Create gauss(x,f(y),s) - RooGaussian model_3("model_3", "Gaussian with shifting mean", x, fy_3, sigma); - - // U s i n g R o o P r o d u c t t o t a i l o r p d f - // ----------------------------------------------------------------------- - - // Create product function f(y) = a1*y - RooProduct fy_4("fy_4", "a1*y", RooArgSet(a1, y)); - - // Create gauss(x,f(y),s) - RooGaussian model_4("model_4", "Gaussian with shifting mean", x, fy_4, sigma); - - // P l o t a l l p d f s - // ---------------------------- - - // Make two-dimensional plots in x vs y - TH1 *hh_model_1 = model_1.createHistogram("hh_model_1", x, Binning(50), YVar(y, Binning(50))); - TH1 *hh_model_2 = model_2.createHistogram("hh_model_2", x, Binning(50), YVar(y, Binning(50))); - TH1 *hh_model_3 = model_3.createHistogram("hh_model_3", x, Binning(50), YVar(y, Binning(50))); - TH1 *hh_model_4 = model_4.createHistogram("hh_model_4", x, Binning(50), YVar(y, Binning(50))); - hh_model_1->SetLineColor(kBlue); - hh_model_2->SetLineColor(kBlue); - hh_model_3->SetLineColor(kBlue); - hh_model_4->SetLineColor(kBlue); - - regTH(hh_model_1, "rf202_model2d_1"); - regTH(hh_model_2, "rf202_model2d_2"); - regTH(hh_model_3, "rf202_model2d_3"); - regTH(hh_model_4, "rf202_model2d_4"); - - return true; - } -}; - // Use of tailored p.d.f as conditional p.d.fs.s // // pdf = gauss(x,f(y),sx | y ) with f(y) = a0 + a1*y @@ -1821,86 +1519,6 @@ class TestBasic307 : public RooUnitTest { } }; -// Examples on normalization of p.d.f.s, integration of p.d.fs, construction of -// cumulative distribution functions from p.d.f.s in two dimensions. -class TestBasic308 : public RooUnitTest { -public: - TestBasic308(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Normalization of p.d.f.s in 2D", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // S e t u p m o d e l - // --------------------- - - // Create observables x,y - RooRealVar x("x", "x", -10, 10); - RooRealVar y("y", "y", -10, 10); - - // Create p.d.f. gaussx(x,-2,3), gaussy(y,2,2) - RooGaussian gx("gx", "gx", x, -2.0, 3.0); - RooGaussian gy("gy", "gy", y, +2.0, 2.0); - - // Create gxy = gx(x)*gy(y) - RooProdPdf gxy("gxy", "gxy", RooArgSet(gx, gy)); - - // R e t r i e v e r a w & n o r m a l i z e d v a l u e s o f R o o F i t p . d . f . s - // -------------------------------------------------------------------------------------------------- - - // Return 'raw' unnormalized value of gx - regValue(gxy.getVal(), "rf308_gxy"); - - // Return value of gxy normalized over x _and_ y in range [-10,10] - RooArgSet nset_xy(x, y); - regValue(gxy.getVal(&nset_xy), "rf308_gx_Norm[x,y]"); - - // Create object representing integral over gx - // which is used to calculate gx_Norm[x,y] == gx / gx_Int[x,y] - std::unique_ptr igxy{gxy.createIntegral(RooArgSet(x, y))}; - regValue(igxy->getVal(), "rf308_gx_Int[x,y]"); - - // NB: it is also possible to do the following - - // Return value of gxy normalized over x in range [-10,10] (i.e. treating y as parameter) - RooArgSet nset_x(x); - regValue(gxy.getVal(&nset_x), "rf308_gx_Norm[x]"); - - // Return value of gxy normalized over y in range [-10,10] (i.e. treating x as parameter) - RooArgSet nset_y(y); - regValue(gxy.getVal(&nset_y), "rf308_gx_Norm[y]"); - - // I n t e g r a t e n o r m a l i z e d p d f o v e r s u b r a n g e - // ---------------------------------------------------------------------------- - - // Define a range named "signal" in x from -5,5 - x.setRange("signal", -5, 5); - y.setRange("signal", -3, 3); - - // Create an integral of gxy_Norm[x,y] over x and y in range "signal" - // This is the fraction of of p.d.f. gxy_Norm[x,y] which is in the - // range named "signal" - std::unique_ptr igxy_sig{ - gxy.createIntegral(RooArgSet(x, y), NormSet(RooArgSet(x, y)), Range("signal"))}; - regValue(igxy_sig->getVal(), "rf308_gx_Int[x,y|signal]_Norm[x,y]"); - - // C o n s t r u c t c u m u l a t i v e d i s t r i b u t i o n f u n c t i o n f r o m p d f - // ----------------------------------------------------------------------------------------------------- - - // Create the cumulative distribution function of gx - // i.e. calculate Int[-10,x] gx(x') dx' - std::unique_ptr gxy_cdf{gxy.createCdf(RooArgSet(x, y))}; - - // Plot cdf of gx versus x - TH1 *hh_cdf = gxy_cdf->createHistogram("hh_cdf", x, Binning(40), YVar(y, Binning(40))); - - regTH(hh_cdf, "rf308_cdf"); - - return true; - } -}; - // Projecting p.d.f and data slices in discrete observables. class TestBasic310 : public RooUnitTest { public: @@ -2206,70 +1824,6 @@ class TestBasic313 : public RooUnitTest { } }; -// Working with parameterized ranges in a fit. This an example of a -// fit with an acceptance that changes per-event -// -// pdf = exp(-t/tau) with t[tmin,5] -// -// where t and tmin are both observables in the dataset -class TestBasic314 : public RooUnitTest { -public: - TestBasic314(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Fit with non-rectangular observable boundaries", refFile, writeRef, verbose) - { - } - bool isTestAvailable() override { return !useCodegenBackend(); } - bool testCode() override - { - - // D e f i n e o b s e r v a b l e s a n d d e c a y p d f - // --------------------------------------------------------------- - - // Declare observables - RooRealVar t("t", "t", 0, 5); - RooRealVar tmin("tmin", "tmin", 0, 0, 5); - - // Make parameterized range in t : [tmin,5] - t.setRange(tmin, RooConst(t.getMax())); - - // Make pdf - RooRealVar tau("tau", "tau", -1.54, -10, -0.1); - RooExponential model("model", "model", t, tau); - - // C r e a t e i n p u t d a t a - // ------------------------------------ - - // Generate complete dataset without acceptance cuts (for reference) - std::unique_ptr dall{model.generate(t, 10000)}; - - // Generate a (fake) prototype dataset for acceptance limit values - std::unique_ptr tmp{RooGaussian("gmin", "gmin", tmin, 0.0, 0.5).generate(tmin, 5000)}; - - // Generate dataset with t values that observe (t>tmin) - std::unique_ptr dacc{model.generate(t, ProtoData(*tmp))}; - - // F i t p d f t o d a t a i n a c c e p t a n c e r e g i o n - // ----------------------------------------------------------------------- - - std::unique_ptr r{model.fitTo(*dacc, Save())}; - - // P l o t f i t t e d p d f o n f u l l a n d a c c e p t e d d a t a - // --------------------------------------------------------------------------------- - - // Make plot frame, add datasets and overlay model - RooPlot *frame = t.frame(Title("Fit to data with per-event acceptance")); - dall->plotOn(frame, MarkerColor(kRed), LineColor(kRed)); - model.plotOn(frame); - dacc->plotOn(frame, Name("dacc")); - - // Print fit results to demonstrate absence of bias - regResult(std::move(r), "rf314_fit"); - regPlot(frame, "rf314_plot1"); - - return true; - } -}; - // Marginizalization of multi-dimensional p.d.f.s through integration. class TestBasic315 : public RooUnitTest { public: @@ -2426,100 +1980,6 @@ class TestBasic316 : public RooUnitTest { } }; -// Tools for manipulation of (un)binned datasets -class TestBasic402 : public RooUnitTest { -public: - TestBasic402(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Basic operations on datasets", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // Binned (RooDataHist) and unbinned datasets (RooDataSet) share - // many properties and inherit from a common abstract base class - // (RooAbsData), that provides an interface for all operations - // that can be performed regardless of the data format - - RooRealVar x("x", "x", -10, 10); - RooRealVar y("y", "y", 0, 40); - RooCategory c("c", "c"); - c.defineType("Plus", +1); - c.defineType("Minus", -1); - - // B a s i c O p e r a t i o n s o n u n b i n n e d d a t a s e t s - // -------------------------------------------------------------- - - // RooDataSet is an unbinned dataset (a collection of points in N-dimensional space) - RooDataSet d("d", "d", RooArgSet(x, y, c)); - - // Unlike RooAbsArgs (RooAbsPdf,RooFormulaVar,....) datasets are not attached to - // the variables they are constructed from. Instead they are attached to an internal - // clone of the supplied set of arguments - - // Fill d with dummy values - int i; - for (i = 0; i < 1000; i++) { - x = i / 50 - 10; - y = sqrt(1.0 * i); - c.setLabel((i % 2) ? "Plus" : "Minus"); - - // We must explicitly refer to x,y,c here to pass the values because - // d is not linked to them (as explained above) - d.add(RooArgSet(x, y, c)); - } - - // R e d u c i n g , A p p e n d i n g a n d M e r g i n g - // ------------------------------------------------------------- - - // The reduce() function returns a new dataset which is a subset of the original - std::unique_ptr d1{d.reduce(SelectVars({x, c}))}; - std::unique_ptr d2{d.reduce(SelectVars(y))}; - std::unique_ptr d3{d.reduce(Cut("y>5.17"))}; - std::unique_ptr d4{d.reduce(SelectVars({x, c}), Cut("y>5.17"))}; - - regValue(d3->numEntries(), "rf403_nd3"); - regValue(d4->numEntries(), "rf403_nd4"); - - // The merge() function adds two data set column-wise - static_cast(*d1).merge(static_cast(d2.get())); - - // The append() function adds two datasets row-wise - static_cast(*d1).append(static_cast(*d3)); - - regValue(d1->numEntries(), "rf403_nd1"); - - // O p e r a t i o n s o n b i n n e d d a t a s e t s - // --------------------------------------------------------- - - // A binned dataset can be constructed empty, from an unbinned dataset, or - // from a ROOT native histogram (TH1,2,3) - - // The binning of real variables (like x,y) is done using their fit range - // 'get/setRange()' and number of specified fit bins 'get/setBins()'. - // Category dimensions of binned datasets get one bin per defined category state - x.setBins(10); - y.setBins(10); - RooDataHist dh("dh", "binned version of d", RooArgSet(x, y), d); - - RooPlot *yframe = y.frame(Bins(10), Title("Operations on binned datasets")); - dh.plotOn(yframe); // plot projection of 2D binned data on y - - // Reduce the 2-dimensional binned dataset to a 1-dimensional binned dataset - // - // All reduce() methods are interfaced in RooAbsData. All reduction techniques - // demonstrated on unbinned datasets can be applied to binned datasets as well. - std::unique_ptr dh2{dh.reduce(SelectVars(y), Cut("x>0"))}; - - // Add dh2 to yframe and redraw - dh2->plotOn(yframe, LineColor(kRed), MarkerColor(kRed), Name("dh2")); - - regPlot(yframe, "rf402_plot1"); - - return true; - } -}; - // Using weights in unbinned datasets. class TestBasic403 : public RooUnitTest { public: @@ -2639,250 +2099,6 @@ class TestBasic403 : public RooUnitTest { } }; -// Working with RooCategory objects to describe discrete variables. -class TestBasic404 : public RooUnitTest { -public: - TestBasic404(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Categories basic functionality", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // C o n s t r u c t a c a t e g o r y w i t h l a b e l s - // ---------------------------------------------------------------- - - // Define a category with labels only - RooCategory tagCat("tagCat", "Tagging category"); - tagCat.defineType("Lepton"); - tagCat.defineType("Kaon"); - tagCat.defineType("NetTagger-1"); - tagCat.defineType("NetTagger-2"); - - // C o n s t r u c t a c a t e g o r y w i t h l a b e l s a n d i n d e c e s - // ---------------------------------------------------------------------------------------- - - // Define a category with explicitly numbered states - RooCategory b0flav("b0flav", "B0 flavour eigenstate"); - b0flav.defineType("B0", -1); - b0flav.defineType("B0bar", 1); - - // G e n e r a t e d u m m y d a t a f o r t a b u l a t i o n d e m o - // ---------------------------------------------------------------------------- - - // Generate a dummy dataset - RooRealVar x("x", "x", 0, 10); - std::unique_ptr data{RooPolynomial("p", "p", x).generate({x, b0flav, tagCat}, 10000)}; - - // P r i n t t a b l e s o f c a t e g o r y c o n t e n t s o f d a t a s e t s - // ------------------------------------------------------------------------------------------ - - // Tables are equivalent of plots for categories - Roo1DTable *btable = data->table(b0flav); - - // Create table for subset of events matching cut expression - Roo1DTable *ttable = data->table(tagCat, "x>8.23"); - - // Create table for all (tagCat x b0flav) state combinations - Roo1DTable *bttable = data->table(RooArgSet(tagCat, b0flav)); - - // Retrieve number of events from table - // Number can be non-integer if source dataset has weighed events - double nb0 = btable->get("B0"); - regValue(nb0, "rf404_nb0"); - - // Retrieve fraction of events with "Lepton" tag - double fracLep = ttable->getFrac("Lepton"); - regValue(fracLep, "rf404_fracLep"); - - // D e f i n i n g r a n g e s f o r p l o t t i n g , f i t t i n g o n c a t e g o r i e s - // ------------------------------------------------------------------------------------------------------ - - // Define named range as comma separated list of labels - tagCat.setRange("good", "Lepton,Kaon"); - - // Or add state names one by one - tagCat.addToRange("soso", "NetTagger-1"); - tagCat.addToRange("soso", "NetTagger-2"); - - // Use category range in dataset reduction specification - std::unique_ptr goodData{data->reduce(CutRange("good"))}; - Roo1DTable *gtable = goodData->table(tagCat); - - regTable(btable, "rf404_btable"); - regTable(ttable, "rf404_ttable"); - regTable(bttable, "rf404_bttable"); - regTable(gtable, "rf404_gtable"); - - return true; - } -}; - -// Demonstration of real-->discrete mapping functions. -class TestBasic405 : public RooUnitTest { -public: - TestBasic405(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Real-to-category functions", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // D e f i n e p d f i n x , s a m p l e d a t a s e t i n x - // ------------------------------------------------------------------------ - - // Define a dummy PDF in x - RooRealVar x("x", "x", 0, 10); - RooArgusBG a("a", "argus(x)", x, RooRealConstant::value(10), RooRealConstant::value(-1)); - - // Generate a dummy dataset - std::unique_ptr data{a.generate(x, 10000)}; - - // C r e a t e a t h r e s h o l d r e a l - > c a t f u n c t i o n - // -------------------------------------------------------------------------- - - // A RooThresholdCategory is a category function that maps regions in a real-valued - // input observable observables to state names. At construction time a 'default' - // state name must be specified to which all values of x are mapped that are not - // otherwise assigned - RooThresholdCategory xRegion("xRegion", "region of x", x, "Background"); - - // Specify thresholds and state assignments one-by-one. - // Each statement specifies that all values _below_ the given value - // (and above any lower specified threshold) are mapped to the - // category state with the given name - // - // Background | SideBand | Signal | SideBand | Background - // 4.23 5.23 8.23 9.23 - xRegion.addThreshold(4.23, "Background"); - xRegion.addThreshold(5.23, "SideBand"); - xRegion.addThreshold(8.23, "Signal"); - xRegion.addThreshold(9.23, "SideBand"); - - // U s e t h r e s h o l d f u n c t i o n t o p l o t d a t a r e g i o n s - // ------------------------------------------------------------------------------------- - - // Add values of threshold function to dataset so that it can be used as observable - data->addColumn(xRegion); - - // Make plot of data in x - RooPlot *xframe = x.frame(Title("Demo of threshold and binning mapping functions")); - data->plotOn(xframe); - - // Use calculated category to select sideband data - data->plotOn(xframe, Cut("xRegion==xRegion::SideBand"), MarkerColor(kRed), LineColor(kRed), Name("data_cut")); - - // C r e a t e a b i n n i n g r e a l - > c a t f u n c t i o n - // ---------------------------------------------------------------------- - - // A RooBinningCategory is a category function that maps bins of a (named) binning definition - // in a real-valued input observable observables to state names. The state names are automatically - // constructed from the variable name, the binning name and the bin number. If no binning name - // is specified the default binning is mapped - - x.setBins(10, "coarse"); - RooBinningCategory xBins("xBins", "coarse bins in x", x, "coarse"); - - // U s e b i n n i n g f u n c t i o n f o r t a b u l a t i o n a n d p l o t t i n g - // ----------------------------------------------------------------------------------------------- - - // Print table of xBins state multiplicity. Note that xBins does not need to be an observable in data - // it can be a function of observables in data as well - Roo1DTable *xbtable = data->table(xBins); - - // Add values of xBins function to dataset so that it can be used as observable - RooCategory *xb = (RooCategory *)data->addColumn(xBins); - - // Define range "alt" as including bins 1,3,5,7,9 - xb->setRange("alt", "x_coarse_bin1,x_coarse_bin3,x_coarse_bin5,x_coarse_bin7,x_coarse_bin9"); - - // Construct subset of data matching range "alt" but only for the first 5000 events and plot it on the frame - std::unique_ptr dataSel{data->reduce(CutRange("alt"), EventRange(0, 5000))}; - // dataSel->plotOn(xframe,MarkerColor(kGreen),LineColor(kGreen),Name("data_sel")) ; - - regTable(xbtable, "rf405_xbtable"); - regPlot(xframe, "rf405_plot1"); - - return true; - } -}; - -// Demonstration of discrete-->discrete (invertable) functions. -class TestBasic406 : public RooUnitTest { -public: - TestBasic406(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Category-to-category functions", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // C o n s t r u c t t w o c a t e g o r i e s - // ---------------------------------------------- - - // Define a category with labels only - RooCategory tagCat("tagCat", "Tagging category"); - tagCat.defineType("Lepton"); - tagCat.defineType("Kaon"); - tagCat.defineType("NetTagger-1"); - tagCat.defineType("NetTagger-2"); - - // Define a category with explicitly numbered states - RooCategory b0flav("b0flav", "B0 flavour eigenstate"); - b0flav.defineType("B0", -1); - b0flav.defineType("B0bar", 1); - - // Construct a dummy dataset with random values of tagCat and b0flav - RooRealVar x("x", "x", 0, 10); - RooPolynomial p("p", "p", x); - std::unique_ptr data{p.generate({x, b0flav, tagCat}, 10000)}; - - // C r e a t e a c a t - > c a t m a p p i n g c a t e g o r y - // --------------------------------------------------------------------- - - // A RooMappedCategory is category->category mapping function based on string expression - // The constructor takes an input category an a default state name to which unassigned - // states are mapped - RooMappedCategory tcatType("tcatType", "tagCat type", tagCat, "Cut based"); - - // Enter fully specified state mappings - tcatType.map("Lepton", "Cut based"); - tcatType.map("Kaon", "Cut based"); - - // Enter a wildcard expression mapping - tcatType.map("NetTagger*", "Neural Network"); - - // Make a table of the mapped category state multiplicit in data - Roo1DTable *mtable = data->table(tcatType); - - // C r e a t e a c a t X c a t p r o d u c t c a t e g o r y - // ---------------------------------------------------------------------- - - // A SUPER-category is 'product' of _lvalue_ categories. The state names of a super - // category is a composite of the state labels of the input categories - RooSuperCategory b0Xtcat("b0Xtcat", "b0flav X tagCat", RooArgSet(b0flav, tagCat)); - - // Make a table of the product category state multiplicity in data - Roo1DTable *stable = data->table(b0Xtcat); - - // Since the super category is an lvalue, assignment is explicitly possible - b0Xtcat.setLabel("{B0bar;Lepton}"); - - // A MULTI-category is a 'product' of any category (function). The state names of a super - // category is a composite of the state labels of the input categories - RooMultiCategory b0Xttype("b0Xttype", "b0flav X tagType", RooArgSet(b0flav, tcatType)); - - // Make a table of the product category state multiplicity in data - Roo1DTable *xtable = data->table(b0Xttype); - - regTable(mtable, "rf406_mtable"); - regTable(stable, "rf406_stable"); - regTable(xtable, "rf406_xtable"); - - return true; - } -}; - // Using simultaneous p.d.f.s to describe simultaneous fits to multiple // datasets. class TestBasic501 : public RooUnitTest { @@ -3338,90 +2554,6 @@ class TestBasic604 : public RooUnitTest { } }; -// Working with the profile likelihood estimator. -class TestBasic605 : public RooUnitTest { -public: - TestBasic605(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Profile Likelihood operator", refFile, writeRef, verbose) - { - } - bool isTestAvailable() override { return !useCodegenBackend(); } - bool testCode() override - { - - // C r e a t e m o d e l a n d d a t a s e t - // ----------------------------------------------- - - // Observable - RooRealVar x("x", "x", -20, 20); - - // Model (intentional strong correlations) - RooRealVar mean("mean", "mean of g1 and g2", 0, -10, 10); - RooRealVar sigma_g1("sigma_g1", "width of g1", 3); - RooGaussian g1("g1", "g1", x, mean, sigma_g1); - - RooRealVar sigma_g2("sigma_g2", "width of g2", 4, 3.0, 6.0); - RooGaussian g2("g2", "g2", x, mean, sigma_g2); - - RooRealVar frac("frac", "frac", 0.5, 0.0, 1.0); - RooAddPdf model("model", "model", RooArgList(g1, g2), frac); - - // Generate 1000 events - std::unique_ptr data{model.generate(x, 1000)}; - - // C o n s t r u c t p l a i n l i k e l i h o o d - // --------------------------------------------------- - - // Construct unbinned likelihood - std::unique_ptr nll{model.createNLL(*data)}; - nll->SetName("nll"); - - // Minimize likelihood w.r.t all parameters before making plots - RooMinimizer(*nll).migrad(); - - // Plot likelihood scan frac - RooPlot *frame1 = frac.frame(Bins(10), Range(0.01, 0.95), Title("LL and profileLL in frac")); - nll->plotOn(frame1, ShiftToZero()); - - // Plot likelihood scan in sigma_g2 - RooPlot *frame2 = sigma_g2.frame(Bins(10), Range(3.3, 5.0), Title("LL and profileLL in sigma_g2")); - nll->plotOn(frame2, ShiftToZero()); - - // C o n s t r u c t p r o f i l e l i k e l i h o o d i n f r a c - // ----------------------------------------------------------------------- - - // The profile likelihood estimator on nll for frac will minimize nll w.r.t - // all floating parameters except frac for each evaluation - RooProfileLL pll_frac("pll_frac", "pll_frac", *nll, frac); - - // Plot the profile likelihood in frac - pll_frac.plotOn(frame1, LineColor(kRed)); - - // Adjust frame maximum for visual clarity - frame1->SetMinimum(0); - frame1->SetMaximum(3); - - // C o n s t r u c t p r o f i l e l i k e l i h o o d i n s i g m a _ g 2 - // ------------------------------------------------------------------------------- - - // The profile likelihood estimator on nll for sigma_g2 will minimize nll - // w.r.t all floating parameters except sigma_g2 for each evaluation - RooProfileLL pll_sigmag2("pll_sigmag2", "pll_sigmag2", *nll, sigma_g2); - - // Plot the profile likelihood in sigma_g2 - pll_sigmag2.plotOn(frame2, LineColor(kRed)); - - // Adjust frame maximum for visual clarity - frame2->SetMinimum(0); - frame2->SetMaximum(3); - - regPlot(frame1, "rf605_plot1"); - regPlot(frame2, "rf605_plot2"); - - return true; - } -}; - // Understanding and customizing error handling in likelihood evaluations. class TestBasic606 : public RooUnitTest { public: @@ -3536,312 +2668,6 @@ class TestBasic607 : public RooUnitTest { } }; -// Unbinned maximum likelihood fit of an efficiency eff(x) function to -// a dataset D(x,cut), where cut is a category encoding a selection, of which -// the efficiency as function of x should be described by eff(x). -class TestBasic701 : public RooUnitTest { -public: - TestBasic701(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Efficiency operator p.d.f. 1D", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - // C o n s t r u c t e f f i c i e n c y f u n c t i o n e ( x ) - // ------------------------------------------------------------------- - - // Declare variables x,mean,sigma with associated name, title, initial value and allowed range - RooRealVar x("x", "x", -10, 10); - - // Efficiency function eff(x;a,b) - RooRealVar a("a", "a", 0.4, 0, 1); - RooRealVar b("b", "b", 5); - RooRealVar c("c", "c", -1, -10, 10); - RooFormulaVar effFunc("effFunc", "(1-a)+a*cos((x-c)/b)", RooArgList(a, b, c, x)); - - // C o n s t r u c t c o n d i t i o n a l e f f i c i e n c y p d f E ( c u t | x ) - // ------------------------------------------------------------------------------------------ - - // Acceptance state cut (1 or 0) - RooCategory cut("cut", "cutr"); - cut.defineType("accept", 1); - cut.defineType("reject", 0); - - // Construct efficiency p.d.f eff(cut|x) - RooEfficiency effPdf("effPdf", "effPdf", effFunc, cut, "accept"); - - // G e n e r a t e d a t a ( x , c u t ) f r o m a t o y m o d e l - // ----------------------------------------------------------------------------- - - // Construct global shape p.d.f shape(x) and product model(x,cut) = eff(cut|x)*shape(x) - // (These are _only_ needed to generate some toy MC here to be used later) - RooPolynomial shapePdf("shapePdf", "shapePdf", x, RooConst(-0.095)); - RooProdPdf model("model", "model", shapePdf, Conditional(effPdf, cut)); - - // Generate some toy data from model - std::unique_ptr data{model.generate({x, cut}, 10000)}; - - // F i t c o n d i t i o n a l e f f i c i e n c y p d f t o d a t a - // -------------------------------------------------------------------------- - - // Fit conditional efficiency p.d.f to data - effPdf.fitTo(*data, ConditionalObservables(x)); - - // P l o t f i t t e d , d a t a e f f i c i e n c y - // -------------------------------------------------------- - - // Plot distribution of all events and accepted fraction of events on frame - RooPlot *frame1 = x.frame(Bins(20), Title("Data (all, accepted)")); - data->plotOn(frame1); - data->plotOn(frame1, Cut("cut==cut::accept"), MarkerColor(kRed), LineColor(kRed)); - - // Plot accept/reject efficiency on data overlay fitted efficiency curve - RooPlot *frame2 = x.frame(Bins(20), Title("Fitted efficiency")); - data->plotOn(frame2, Efficiency(cut)); // needs ROOT version >= 5.21 - effFunc.plotOn(frame2, LineColor(kRed)); - - regPlot(frame1, "rf701_plot1"); - regPlot(frame2, "rf701_plot2"); - - return true; - } -}; - -// Unbinned maximum likelihood fit of an efficiency eff(x) function to -// a dataset D(x,cut), where cut is a category encoding a selection whose -// efficiency as function of x should be described by eff(x). -class TestBasic702 : public RooUnitTest { -public: - TestBasic702(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Efficiency operator p.d.f. 2D", refFile, writeRef, verbose) - { - } - bool testCode() override - { - - bool flat = false; - - // C o n s t r u c t e f f i c i e n c y f u n c t i o n e ( x , y ) - // ----------------------------------------------------------------------- - - // Declare variables x,mean,sigma with associated name, title, initial value and allowed range - RooRealVar x("x", "x", -10, 10); - RooRealVar y("y", "y", -10, 10); - - // Efficiency function eff(x;a,b) - RooRealVar ax("ax", "ay", 0.6, 0, 1); - RooRealVar bx("bx", "by", 5); - RooRealVar cx("cx", "cy", -1, -10, 10); - - RooRealVar ay("ay", "ay", 0.2, 0, 1); - RooRealVar by("by", "by", 5); - RooRealVar cy("cy", "cy", -1, -10, 10); - - RooFormulaVar effFunc("effFunc", "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))", - RooArgList(ax, bx, cx, x, ay, by, cy, y)); - - // Acceptance state cut (1 or 0) - RooCategory cut("cut", "cutr"); - cut.defineType("accept", 1); - cut.defineType("reject", 0); - - // C o n s t r u c t c o n d i t i o n a l e f f i c i e n c y p d f E ( c u t | x , y ) - // --------------------------------------------------------------------------------------------- - - // Construct efficiency p.d.f eff(cut|x) - RooEfficiency effPdf("effPdf", "effPdf", effFunc, cut, "accept"); - - // G e n e r a t e d a t a ( x , y , c u t ) f r o m a t o y m o d e l - // ------------------------------------------------------------------------------- - - // Construct global shape p.d.f shape(x) and product model(x,cut) = eff(cut|x)*shape(x) - // (These are _only_ needed to generate some toy MC here to be used later) - RooPolynomial shapePdfX("shapePdfX", "shapePdfX", x, RooConst(flat ? 0 : -0.095)); - RooPolynomial shapePdfY("shapePdfY", "shapePdfY", y, RooConst(flat ? 0 : +0.095)); - RooProdPdf shapePdf("shapePdf", "shapePdf", RooArgSet(shapePdfX, shapePdfY)); - RooProdPdf model("model", "model", shapePdf, Conditional(effPdf, cut)); - - // Generate some toy data from model - std::unique_ptr data{model.generate({x, y, cut}, 10000)}; - - // F i t c o n d i t i o n a l e f f i c i e n c y p d f t o d a t a - // -------------------------------------------------------------------------- - - // Fit conditional efficiency p.d.f to data - effPdf.fitTo(*data, ConditionalObservables(RooArgSet(x, y))); - - // P l o t f i t t e d , d a t a e f f i c i e n c y - // -------------------------------------------------------- - - // Make 2D histograms of all data, selected data and efficiency function - TH1 *hh_data_all = data->createHistogram("hh_data_all", x, Binning(8), YVar(y, Binning(8))); - TH1 *hh_data_sel = - data->createHistogram("hh_data_sel", x, Binning(8), YVar(y, Binning(8)), Cut("cut==cut::accept")); - TH1 *hh_eff = effFunc.createHistogram("hh_eff", x, Binning(50), YVar(y, Binning(50))); - - // Some adjustsment for good visualization - hh_data_all->SetMinimum(0); - hh_data_sel->SetMinimum(0); - hh_eff->SetMinimum(0); - hh_eff->SetLineColor(kBlue); - - regTH(hh_data_all, "rf702_hh_data_all"); - regTH(hh_data_sel, "rf702_hh_data_sel"); - regTH(hh_eff, "rf702_hh_eff"); - - return true; - } -}; - -// Using a product of an (acceptance) efficiency and a p.d.f as p.d.f. -class TestBasic703 : public RooUnitTest { -public: - TestBasic703(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Efficiency product operator p.d.f", refFile, writeRef, verbose) - { - } - bool isTestAvailable() override { return !useCodegenBackend(); } - bool testCode() override - { - - // D e f i n e o b s e r v a b l e s a n d d e c a y p d f - // --------------------------------------------------------------- - - // Declare observables - RooRealVar t("t", "t", 0, 5); - - // Make pdf - RooRealVar tau("tau", "tau", -1.54, -4, -0.1); - RooExponential model("model", "model", t, tau); - - // D e f i n e e f f i c i e n c y f u n c t i o n - // --------------------------------------------------- - - // Use error function to simulate turn-on slope - RooFormulaVar eff("eff", "0.5*(std::erf((t-1)/0.5)+1)", t); - - // D e f i n e d e c a y p d f w i t h e f f i c i e n c y - // --------------------------------------------------------------- - - // Multiply pdf(t) with efficiency in t - RooEffProd modelEff("modelEff", "model with efficiency", model, eff); - - // P l o t e f f i c i e n c y , p d f - // ---------------------------------------- - - RooPlot *frame1 = t.frame(Title("Efficiency")); - eff.plotOn(frame1, LineColor(kRed)); - - RooPlot *frame2 = t.frame(Title("Pdf with and without efficiency")); - - model.plotOn(frame2, LineStyle(kDashed)); - modelEff.plotOn(frame2); - - // G e n e r a t e t o y d a t a , f i t m o d e l E f f t o d a t a - // ------------------------------------------------------------------------------ - - // Generate events. If the input pdf has an internal generator, the internal generator - // is used and an accept/reject sampling on the efficiency is applied. - std::unique_ptr data{modelEff.generate(t, 10000)}; - - // Fit pdf. The normalization integral is calculated numerically. - modelEff.fitTo(*data); - - // Plot generated data and overlay fitted pdf - RooPlot *frame3 = t.frame(Title("Fitted pdf with efficiency")); - data->plotOn(frame3); - modelEff.plotOn(frame3); - - regPlot(frame1, "rf703_plot1"); - regPlot(frame2, "rf703_plot2"); - regPlot(frame3, "rf703_plot3"); - - return true; - } -}; - -// Using a p.d.f defined by a sum of real-valued amplitude components. -class TestBasic704 : public RooUnitTest { -public: - TestBasic704(TFile *refFile, bool writeRef, int verbose) - : RooUnitTest("Amplitude sum operator p.d.f", refFile, writeRef, verbose) - { - } - bool isTestAvailable() override { return !useCodegenBackend(); } - bool testCode() override - { - - // S e t u p 2 D a m p l i t u d e f u n c t i o n s - // ------------------------------------------------------- - - // Observables - RooRealVar t("t", "time", -1., 15.); - RooRealVar cosa("cosa", "cos(alpha)", -1., 1.); - - RooRealVar tau("tau", "#tau", 1.5); - RooRealVar deltaGamma("deltaGamma", "deltaGamma", 0.3); - RooFormulaVar coshG("coshGBasis", "exp(-@0/ @1)*cosh(@0*@2/2)", {t, tau, deltaGamma}); - RooFormulaVar sinhG("sinhGBasis", "exp(-@0/ @1)*sinh(@0*@2/2)", {t, tau, deltaGamma}); - - // Construct polynomial amplitudes in cos(a) - RooPolyVar poly1("poly1", "poly1", cosa, RooArgList(0.5, 0.2, 0.2), 0); - RooPolyVar poly2("poly2", "poly2", cosa, RooArgList(1.0, -0.2, 3.0), 0); - - // Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa) - RooProduct ampl1("ampl1", "amplitude 1", {poly1, coshG}); - RooProduct ampl2("ampl2", "amplitude 2", {poly2, sinhG}); - - // C o n s t r u c t a m p l i t u d e s u m p d f - // ----------------------------------------------------- - - // Amplitude strengths - RooRealVar f1("f1", "f1", 1, 0, 2); - RooRealVar f2("f2", "f2", 0.5, 0, 2); - - // Construct pdf - RooRealSumPdf pdf("pdf", "pdf", RooArgList(ampl1, ampl2), RooArgList(f1, f2)); - - // Generate some toy data from pdf - std::unique_ptr data{pdf.generate({t, cosa}, 10000)}; - - // Fit pdf to toy data with only amplitude strength floating - pdf.fitTo(*data); - - // P l o t a m p l i t u d e s u m p d f - // ------------------------------------------- - - // Make 2D plots of amplitudes - TH1 *hh_cos = ampl1.createHistogram("hh_cos", t, Binning(50), YVar(cosa, Binning(50))); - TH1 *hh_sin = ampl2.createHistogram("hh_sin", t, Binning(50), YVar(cosa, Binning(50))); - hh_cos->SetLineColor(kBlue); - hh_sin->SetLineColor(kBlue); - - // Make projection on t, plot data, pdf and its components - // Note component projections may be larger than sum because amplitudes can be negative - RooPlot *frame1 = t.frame(); - data->plotOn(frame1); - pdf.plotOn(frame1); - pdf.plotOn(frame1, Components(ampl1), LineStyle(kDashed)); - pdf.plotOn(frame1, Components(ampl2), LineStyle(kDashed), LineColor(kRed)); - - // Make projection on cosa, plot data, pdf and its components - // Note that components projection may be larger than sum because amplitudes can be negative - RooPlot *frame2 = cosa.frame(); - data->plotOn(frame2); - pdf.plotOn(frame2); - pdf.plotOn(frame2, Components(ampl1), LineStyle(kDashed)); - pdf.plotOn(frame2, Components(ampl2), LineStyle(kDashed), LineColor(kRed)); - - regPlot(frame1, "rf704_plot1"); - regPlot(frame2, "rf704_plot2"); - regTH(hh_cos, "rf704_hh_cos"); - regTH(hh_sin, "rf704_hh_sin"); - - return true; - } -}; - // Histogram based p.d.f.s and functions. class TestBasic706 : public RooUnitTest { public: diff --git a/roofit/roofitcore/test/testRooAbsPdf.cxx b/roofit/roofitcore/test/testRooAbsPdf.cxx index f13d2afa9e316..cd7b57a6c17de 100644 --- a/roofit/roofitcore/test/testRooAbsPdf.cxx +++ b/roofit/roofitcore/test/testRooAbsPdf.cxx @@ -8,11 +8,14 @@ #include #include #include +#include #include #include #include #include #include +#include +#include #include #include #include @@ -23,11 +26,15 @@ #include #include +#include #include #include "gtest_wrapper.h" +#include #include +#include +#include class FitTest : public testing::TestWithParam> { public: @@ -520,6 +527,169 @@ TEST(RooAbsPdf, NormSetChange) EXPECT_NE(v1, v2); } +namespace { + +/// Integral of an unnormalized Gaussian exp(-0.5 ((x - mean) / sigma)^2) over [lo, hi]. +double gaussInt(double lo, double hi, double mean, double sigma) +{ + const double sqrt2 = std::sqrt(2.0); + return sigma * std::sqrt(TMath::Pi() / 2.) * + (std::erf((hi - mean) / (sqrt2 * sigma)) - std::erf((lo - mean) / (sqrt2 * sigma))); +} + +} // namespace + +/// Normalization, integration and cdf of a pdf in one dimension, checked +/// against the analytically known Gaussian integrals. Replaces the former +/// stressRooFit test based on the rf110 tutorial, which compared against +/// stored reference values. +TEST(RooAbsPdf, Normalization1D) +{ + RooRealVar x("x", "x", -10, 10); + RooGaussian gx("gx", "gx", x, -2.0, 3.0); + + const double rawVal = std::exp(-0.5 * std::pow((x.getVal() + 2.) / 3., 2)); + const double normInt = gaussInt(-10, 10, -2., 3.); + + // Raw unnormalized value and value normalized over x in [-10, 10] + EXPECT_NEAR(gx.getVal(), rawVal, 1e-10); + RooArgSet nset{x}; + EXPECT_NEAR(gx.getVal(&nset), rawVal / normInt, 1e-10); + + // Integral over the full range + std::unique_ptr igx{gx.createIntegral(x)}; + EXPECT_NEAR(igx->getVal(), normInt, 1e-6 * normInt); + + // Fraction of the normalized pdf contained in the "signal" sub range + x.setRange("signal", -5, 5); + std::unique_ptr igxSig{gx.createIntegral(x, RooFit::NormSet(x), RooFit::Range("signal"))}; + const double sigFrac = gaussInt(-5, 5, -2., 3.) / normInt; + EXPECT_NEAR(igxSig->getVal(), sigFrac, 1e-6); + + // Cumulative distribution function + std::unique_ptr cdf{gx.createCdf(x)}; + for (double xVal : {-10., -5., -2., 0., 3., 10.}) { + x.setVal(xVal); + EXPECT_NEAR(cdf->getVal(), gaussInt(-10, xVal, -2., 3.) / normInt, 1e-6) << "cdf at x = " << xVal; + } +} + +/// Normalization and integration of a product pdf in two dimensions, checked +/// against the analytically known Gaussian integrals. Replaces the former +/// stressRooFit test based on the rf308 tutorial, which compared against +/// stored reference values. +TEST(RooAbsPdf, Normalization2D) +{ + RooRealVar x("x", "x", -10, 10); + RooRealVar y("y", "y", -10, 10); + + RooGaussian gx("gx", "gx", x, -2.0, 3.0); + RooGaussian gy("gy", "gy", y, +2.0, 2.0); + RooProdPdf gxy("gxy", "gxy", RooArgSet(gx, gy)); + + const double rawX = std::exp(-0.5 * std::pow((x.getVal() + 2.) / 3., 2)); + const double rawY = std::exp(-0.5 * std::pow((y.getVal() - 2.) / 2., 2)); + const double intX = gaussInt(-10, 10, -2., 3.); + const double intY = gaussInt(-10, 10, +2., 2.); + + EXPECT_NEAR(gxy.getVal(), rawX * rawY, 1e-10); + + // Normalized over both, or only one of the observables (the other one is + // then treated as a parameter) + RooArgSet nsetXY{x, y}; + RooArgSet nsetX{x}; + RooArgSet nsetY{y}; + EXPECT_NEAR(gxy.getVal(&nsetXY), rawX * rawY / (intX * intY), 1e-10); + EXPECT_NEAR(gxy.getVal(&nsetX), rawX / intX, 1e-10); + EXPECT_NEAR(gxy.getVal(&nsetY), rawY / intY, 1e-10); + + std::unique_ptr igxy{gxy.createIntegral({x, y})}; + EXPECT_NEAR(igxy->getVal(), intX * intY, 1e-6 * intX * intY); + + // Fraction of the normalized pdf contained in the rectangular "signal" range + x.setRange("signal", -5, 5); + y.setRange("signal", -3, 3); + std::unique_ptr igxySig{gxy.createIntegral({x, y}, RooFit::NormSet(RooArgSet{x, y}), RooFit::Range("signal"))}; + const double sigFrac = gaussInt(-5, 5, -2., 3.) * gaussInt(-3, 3, 2., 2.) / (intX * intY); + EXPECT_NEAR(igxySig->getVal(), sigFrac, 1e-6); + + // The cdf of the product of two independent pdfs factorizes into the + // product of the marginal cdfs + std::unique_ptr cdf{gxy.createCdf({x, y})}; + const std::vector> cdfPoints{{-5., -2.}, {0., 0.}, {2.5, 4.}, {10., 10.}}; + for (auto const &[xVal, yVal] : cdfPoints) { + x.setVal(xVal); + y.setVal(yVal); + const double ref = gaussInt(-10, xVal, -2., 3.) / intX * gaussInt(-10, yVal, 2., 2.) / intY; + EXPECT_NEAR(cdf->getVal(), ref, 1e-6) << "cdf at (x, y) = (" << xVal << ", " << yVal << ")"; + } +} + +/// Configuration of numeric integration, validated against the analytical +/// integral of the Landau pdf instead of stored reference values. Replaces the +/// former stressRooFit test based on the rf111 tutorial. +TEST(RooAbsPdf, NumIntConfig) +{ + RooRealVar x("x", "x", -10, 10); + RooLandau landau("landau", "landau", x, 0.0, 0.1); + + // The analytical integral serves as the reference + const double refVal = std::unique_ptr{landau.createIntegral(x)}->getVal(); + + // Disable analytic integration and integrate with the default numeric + // integrator configuration + landau.forceNumInt(true); + const double val1 = std::unique_ptr{landau.createIntegral(x)}->getVal(); + EXPECT_NEAR(val1, refVal, 1e-3 * refVal); + + // Use a custom configuration, once passed explicitly to createIntegral() + // and once set as the default configuration of the pdf object + RooNumIntConfig customConfig(*RooAbsReal::defaultIntegratorConfig()); + customConfig.setEpsAbs(1e-8); + customConfig.setEpsRel(1e-8); + + const double val2 = + std::unique_ptr{landau.createIntegral(x, RooFit::NumIntConfig(customConfig))}->getVal(); + EXPECT_NEAR(val2, refVal, 1e-3 * refVal); + + landau.setIntegratorConfig(customConfig); + const double val3 = std::unique_ptr{landau.createIntegral(x)}->getVal(); + + // Both ways of passing the custom configuration must give the identical result + EXPECT_DOUBLE_EQ(val3, val2); +} + +/// Unbinned fit with a per-event acceptance region, implemented via a range +/// that is parameterized by another observable in the dataset. The fit must +/// recover the generating decay constant without bias. Replaces the former +/// stressRooFit test based on the rf314 tutorial. +TEST_P(FitTest, ParameterizedRangeFit) +{ + using namespace RooFit; + + RooRealVar t("t", "t", 0, 5); + RooRealVar tmin("tmin", "tmin", 0, 0, 5); + + // Parameterized range in t : [tmin, 5] + t.setRange(tmin, RooConst(t.getMax())); + + RooRealVar tau("tau", "tau", -1.54, -10, -0.1); + RooExponential model("model", "model", t, tau); + + // Prototype dataset with per-event acceptance limit values + RooGaussian gmin("gmin", "gmin", tmin, 0.0, 0.5); + std::unique_ptr proto{gmin.generate(tmin, 5000)}; + + // Dataset with t values that observe t > tmin + std::unique_ptr data{model.generate(t, ProtoData(*proto))}; + + std::unique_ptr res{model.fitTo(*data, Save(), PrintLevel(-1), _evalBackend)}; + + EXPECT_EQ(res->status(), 0); + EXPECT_EQ(res->covQual(), 3); + expectParamNear(*res, "tau", -1.54); +} + INSTANTIATE_TEST_SUITE_P(RooAbsPdf, FitTest, testing::Values(ROOFIT_EVAL_BACKENDS), [](testing::TestParamInfo const ¶mInfo) { std::stringstream ss; diff --git a/roofit/roofitcore/test/testRooAbsReal.cxx b/roofit/roofitcore/test/testRooAbsReal.cxx index b46ad50e34254..083f8583b3d33 100644 --- a/roofit/roofitcore/test/testRooAbsReal.cxx +++ b/roofit/roofitcore/test/testRooAbsReal.cxx @@ -4,21 +4,28 @@ #include #include +#include #include #include #include +#include +#include #include #include +#include +#include #include #include #include #include #include +#include #include #include +#include #include // ROOT-6882: Cannot read from ULong64_t branches. @@ -130,3 +137,62 @@ TEST(RooAbsReal, YieldsHistogram) EXPECT_FLOAT_EQ(c1 * x.getBinWidth(1), c2 * x.getBinWidth(0)) << "relative yield is wrong"; EXPECT_FLOAT_EQ(c1 + c2, totalYield) << "total yield is wrong"; } + +/// The sum and product utility classes (RooFormulaVar, RooPolyVar, +/// RooAddition, RooProduct) used to tailor the parameters of a pdf, checked +/// against directly computed function values. Replaces the former stressRooFit +/// test based on the rf302 tutorial, which compared against stored reference +/// histograms. +TEST(RooAbsReal, UtilityFunctionsToTailorPdfs) +{ + RooRealVar x("x", "x", -5, 5); + RooRealVar y("y", "y", -5, 5); + + RooRealVar a0("a0", "a0", -1.5, -5, 5); + RooRealVar a1("a1", "a1", -0.5, -1, 1); + RooRealVar sigma("sigma", "width of gaussian", 0.5); + + RooFormulaVar fy1("fy_1", "a0-a1*sqrt(10*abs(y))", RooArgSet(y, a0, a1)); + RooPolyVar fy2("fy_2", "fy_2", y, RooArgSet(a0, a1)); + RooAddition fy3("fy_3", "a0+y", RooArgSet(a0, y)); + RooProduct fy4("fy_4", "a1*y", RooArgSet(a1, y)); + + RooGaussian model1("model_1", "Gaussian with shifting mean", x, fy1, sigma); + RooGaussian model2("model_2", "Gaussian with shifting mean", x, fy2, sigma); + RooGaussian model3("model_3", "Gaussian with shifting mean", x, fy3, sigma); + RooGaussian model4("model_4", "Gaussian with shifting mean", x, fy4, sigma); + + // Normalized value of a Gaussian pdf in x with mean "mu", truncated to the + // range of x. + auto gaussPdfVal = [&](double xVal, double mu) { + const double sqrt2 = std::sqrt(2.0); + const double norm = 0.5 * std::sqrt(2. * TMath::Pi()) * sigma.getVal() * + (std::erf((x.getMax() - mu) / (sqrt2 * sigma.getVal())) - + std::erf((x.getMin() - mu) / (sqrt2 * sigma.getVal()))); + return std::exp(-0.5 * std::pow((xVal - mu) / sigma.getVal(), 2)) / norm; + }; + + RooArgSet normSet{x}; + + for (double yVal : {-4.5, -1.2, 0., 0.7, 3.3}) { + y.setVal(yVal); + + const double mu1 = a0.getVal() - a1.getVal() * std::sqrt(10. * std::abs(yVal)); + const double mu2 = a0.getVal() + a1.getVal() * yVal; + const double mu3 = a0.getVal() + yVal; + const double mu4 = a1.getVal() * yVal; + + EXPECT_NEAR(fy1.getVal(), mu1, 1e-12) << "y = " << yVal; + EXPECT_NEAR(fy2.getVal(), mu2, 1e-12) << "y = " << yVal; + EXPECT_NEAR(fy3.getVal(), mu3, 1e-12) << "y = " << yVal; + EXPECT_NEAR(fy4.getVal(), mu4, 1e-12) << "y = " << yVal; + + for (double xVal : {-3., 0., 1.5}) { + x.setVal(xVal); + EXPECT_NEAR(model1.getVal(&normSet), gaussPdfVal(xVal, mu1), 1e-9) << "(x, y) = (" << xVal << ", " << yVal << ")"; + EXPECT_NEAR(model2.getVal(&normSet), gaussPdfVal(xVal, mu2), 1e-9) << "(x, y) = (" << xVal << ", " << yVal << ")"; + EXPECT_NEAR(model3.getVal(&normSet), gaussPdfVal(xVal, mu3), 1e-9) << "(x, y) = (" << xVal << ", " << yVal << ")"; + EXPECT_NEAR(model4.getVal(&normSet), gaussPdfVal(xVal, mu4), 1e-9) << "(x, y) = (" << xVal << ", " << yVal << ")"; + } + } +} diff --git a/roofit/roofitcore/test/testRooCategory.cxx b/roofit/roofitcore/test/testRooCategory.cxx index f2f17e6f6d15e..f96ed2145e6df 100644 --- a/roofit/roofitcore/test/testRooCategory.cxx +++ b/roofit/roofitcore/test/testRooCategory.cxx @@ -1,15 +1,27 @@ // Tests for the RooCategory // Author: Jonas Rembser, CERN 04/2021 +#include +#include #include #include #include #include +#include +#include +#include +#include +#include #include #include +#include +#include +#include +#include + // GitHub issue 10278: RooDataSet incorrectly loads RooCategory values from TTree branch of type Short_t TEST(RooCategory, CategoryDefineMultiState) { @@ -27,3 +39,245 @@ TEST(RooCategory, CategoryDefineMultiState) EXPECT_EQ(static_cast((*data.get(0))["cat"]).getCurrentIndex(), 2); } + +/// Roo1DTable tabulation of category contents of a dataset, and named ranges +/// on categories used in dataset reduction. The tables are checked against +/// counts computed while deterministically filling the dataset. Replaces the +/// former stressRooFit test based on the rf404 tutorial, which compared +/// against stored reference tables. +TEST(RooCategory, TablesAndRanges) +{ + RooCategory tagCat("tagCat", "Tagging category"); + tagCat.defineType("Lepton"); + tagCat.defineType("Kaon"); + tagCat.defineType("NetTagger-1"); + tagCat.defineType("NetTagger-2"); + + RooCategory b0flav("b0flav", "B0 flavour eigenstate"); + b0flav.defineType("B0", -1); + b0flav.defineType("B0bar", 1); + + RooRealVar x("x", "x", 0, 10); + RooDataSet data("data", "data", {x, b0flav, tagCat}); + + const std::vector tagStates{"Lepton", "Kaon", "NetTagger-1", "NetTagger-2"}; + + // Fill the dataset deterministically and keep track of the expected counts + std::map nFlav; + std::map nTagWithCut; + std::map nComb; + std::map nGood; + int nCutTotal = 0; + for (int i = 0; i < 10000; ++i) { + const std::string &tag = tagStates[i % 4]; + const std::string flav = ((i / 2) % 2) ? "B0bar" : "B0"; + const double xVal = (i % 100) / 10.0; + x.setVal(xVal); + tagCat.setLabel(tag.c_str()); + b0flav.setLabel(flav.c_str()); + data.add({x, b0flav, tagCat}); + + ++nFlav[flav]; + ++nComb["{" + tag + ";" + flav + "}"]; + if (xVal > 8.23) { + ++nTagWithCut[tag]; + ++nCutTotal; + } + if (tag == "Lepton" || tag == "Kaon") { + ++nGood[tag]; + } + } + + // Table of a single category + std::unique_ptr btable{data.table(b0flav)}; + EXPECT_DOUBLE_EQ(btable->get("B0"), nFlav["B0"]); + EXPECT_DOUBLE_EQ(btable->get("B0bar"), nFlav["B0bar"]); + + // Table for the subset of events matching a cut expression + std::unique_ptr ttable{data.table(tagCat, "x>8.23")}; + for (const auto &tag : tagStates) { + EXPECT_DOUBLE_EQ(ttable->get(tag.c_str()), nTagWithCut[tag]) << tag; + EXPECT_DOUBLE_EQ(ttable->getFrac(tag.c_str()), double(nTagWithCut[tag]) / nCutTotal) << tag; + } + + // Table of all state combinations of two categories + std::unique_ptr bttable{data.table({tagCat, b0flav})}; + for (const auto &[label, count] : nComb) { + EXPECT_DOUBLE_EQ(bttable->get(label.c_str()), count) << label; + } + + // Named category ranges used in dataset reduction + tagCat.setRange("good", "Lepton,Kaon"); + tagCat.addToRange("soso", "NetTagger-1"); + tagCat.addToRange("soso", "NetTagger-2"); + + std::unique_ptr goodData{data.reduce(RooFit::CutRange("good"))}; + EXPECT_EQ(goodData->numEntries(), nGood["Lepton"] + nGood["Kaon"]); + std::unique_ptr gtable{goodData->table(tagCat)}; + EXPECT_DOUBLE_EQ(gtable->get("Lepton"), nGood["Lepton"]); + EXPECT_DOUBLE_EQ(gtable->get("Kaon"), nGood["Kaon"]); + EXPECT_DOUBLE_EQ(gtable->get("NetTagger-1"), 0.); + EXPECT_DOUBLE_EQ(gtable->get("NetTagger-2"), 0.); + + std::unique_ptr sosoData{data.reduce(RooFit::CutRange("soso"))}; + EXPECT_EQ(sosoData->numEntries(), 10000 - nGood["Lepton"] - nGood["Kaon"]); +} + +/// Real-to-category mapping functions RooThresholdCategory and +/// RooBinningCategory, checked entry-by-entry against the expected mapping. +/// Replaces the former stressRooFit test based on the rf405 tutorial, which +/// compared against stored reference tables and plots. +TEST(RooCategory, RealToCategoryFunctions) +{ + RooRealVar x("x", "x", 0, 10); + + // Threshold mapping: each threshold assigns all values below it (and above + // any lower threshold) to the given state + RooThresholdCategory xRegion("xRegion", "region of x", x, "Background"); + xRegion.addThreshold(4.23, "Background"); + xRegion.addThreshold(5.23, "SideBand"); + xRegion.addThreshold(8.23, "Signal"); + xRegion.addThreshold(9.23, "SideBand"); + + // Binning mapping based on a named binning + x.setBins(10, "coarse"); + RooBinningCategory xBins("xBins", "coarse bins in x", x, "coarse"); + + auto expectedRegion = [](double v) -> std::string { + if (v < 4.23) + return "Background"; + if (v < 5.23) + return "SideBand"; + if (v < 8.23) + return "Signal"; + if (v < 9.23) + return "SideBand"; + return "Background"; + }; + + RooDataSet data("data", "data", x); + + std::map nRegion; + std::map nBin; + int nAltFirstHalf = 0; + const int nEntries = 1000; + for (int i = 0; i < nEntries; ++i) { + // Values chosen such that they never coincide with a threshold or bin + // boundary, where the expected mapping would be floating-point fragile + const double v = 10. * (i + 0.5) / nEntries; + x.setVal(v); + + const std::string region = expectedRegion(v); + const std::string bin = "x_coarse_bin" + std::to_string(int(v)); + EXPECT_STREQ(xRegion.getCurrentLabel(), region.c_str()) << "x = " << v; + EXPECT_STREQ(xBins.getCurrentLabel(), bin.c_str()) << "x = " << v; + + data.add(x); + ++nRegion[region]; + ++nBin[bin]; + if (int(v) % 2 == 1 && i < nEntries / 2) { + ++nAltFirstHalf; + } + } + + // Add the category functions as columns to the dataset and cross-check the + // stored values entry-by-entry + auto *xr = static_cast(data.addColumn(xRegion)); + auto *xb = static_cast(data.addColumn(xBins)); + for (int i = 0; i < nEntries; ++i) { + const RooArgSet *row = data.get(i); + const double v = row->getRealValue("x"); + EXPECT_STREQ(row->getCatLabel("xRegion"), expectedRegion(v).c_str()) << "entry " << i; + EXPECT_STREQ(row->getCatLabel("xBins"), ("x_coarse_bin" + std::to_string(int(v))).c_str()) << "entry " << i; + } + + // Tabulate the computed columns + std::unique_ptr rtable{data.table(*xr)}; + for (const auto &[label, count] : nRegion) { + EXPECT_DOUBLE_EQ(rtable->get(label.c_str()), count) << label; + } + std::unique_ptr btable{data.table(*xb)}; + for (const auto &[label, count] : nBin) { + EXPECT_DOUBLE_EQ(btable->get(label.c_str()), count) << label; + } + + // Use a named range on the computed category column together with an event + // range in a dataset reduction + xb->setRange("alt", "x_coarse_bin1,x_coarse_bin3,x_coarse_bin5,x_coarse_bin7,x_coarse_bin9"); + std::unique_ptr dataSel{data.reduce(RooFit::CutRange("alt"), RooFit::EventRange(0, nEntries / 2))}; + EXPECT_EQ(dataSel->numEntries(), nAltFirstHalf); +} + +/// Category-to-category mapping functions RooMappedCategory, RooSuperCategory +/// and RooMultiCategory, checked state-by-state against the expected mapping. +/// Replaces the former stressRooFit test based on the rf406 tutorial, which +/// compared against stored reference tables. +TEST(RooCategory, CategoryToCategoryFunctions) +{ + RooCategory tagCat("tagCat", "Tagging category"); + tagCat.defineType("Lepton"); + tagCat.defineType("Kaon"); + tagCat.defineType("NetTagger-1"); + tagCat.defineType("NetTagger-2"); + + RooCategory b0flav("b0flav", "B0 flavour eigenstate"); + b0flav.defineType("B0", -1); + b0flav.defineType("B0bar", 1); + + // Category-to-category mapping with explicit and wildcard expressions + RooMappedCategory tcatType("tcatType", "tagCat type", tagCat, "Cut based"); + tcatType.map("Lepton", "Cut based"); + tcatType.map("Kaon", "Cut based"); + tcatType.map("NetTagger*", "Neural Network"); + + auto expectedType = [](const std::string &tag) -> std::string { + return (tag == "Lepton" || tag == "Kaon") ? "Cut based" : "Neural Network"; + }; + + // Product categories of lvalue categories (super) and of any category + // functions (multi) + RooSuperCategory b0Xtcat("b0Xtcat", "b0flav X tagCat", {b0flav, tagCat}); + RooMultiCategory b0Xttype("b0Xttype", "b0flav X tagType", {b0flav, tcatType}); + + const std::vector tagStates{"Lepton", "Kaon", "NetTagger-1", "NetTagger-2"}; + for (const auto &flav : {"B0", "B0bar"}) { + for (const auto &tag : tagStates) { + b0flav.setLabel(flav); + tagCat.setLabel(tag.c_str()); + EXPECT_STREQ(tcatType.getCurrentLabel(), expectedType(tag).c_str()) << tag; + const std::string superLabel = "{" + std::string(flav) + ";" + tag + "}"; + EXPECT_STREQ(b0Xtcat.getCurrentLabel(), superLabel.c_str()); + const std::string multiLabel = "{" + std::string(flav) + ";" + expectedType(tag) + "}"; + EXPECT_STREQ(b0Xttype.getCurrentLabel(), multiLabel.c_str()); + } + } + + // A super category is an lvalue: assigning a state must propagate to the + // input categories + EXPECT_FALSE(b0Xtcat.setLabel("{B0bar;Kaon}")); // returns true on error + EXPECT_STREQ(b0flav.getCurrentLabel(), "B0bar"); + EXPECT_STREQ(tagCat.getCurrentLabel(), "Kaon"); + + // Tabulate the mapped categories in a deterministically filled dataset + RooDataSet data("data", "data", {b0flav, tagCat}); + std::map nType; + std::map nSuper; + for (int i = 0; i < 1000; ++i) { + const std::string &tag = tagStates[i % 4]; + const std::string flav = ((i / 2) % 2) ? "B0bar" : "B0"; + tagCat.setLabel(tag.c_str()); + b0flav.setLabel(flav.c_str()); + data.add({b0flav, tagCat}); + ++nType[expectedType(tag)]; + ++nSuper["{" + flav + ";" + tag + "}"]; + } + + std::unique_ptr mtable{data.table(tcatType)}; + EXPECT_DOUBLE_EQ(mtable->get("Cut based"), nType["Cut based"]); + EXPECT_DOUBLE_EQ(mtable->get("Neural Network"), nType["Neural Network"]); + + std::unique_ptr stable{data.table(b0Xtcat)}; + for (const auto &[label, count] : nSuper) { + EXPECT_DOUBLE_EQ(stable->get(label.c_str()), count) << label; + } +} diff --git a/roofit/roofitcore/test/testRooDataSet.cxx b/roofit/roofitcore/test/testRooDataSet.cxx index 7e1e45443feaa..f4760a752391b 100644 --- a/roofit/roofitcore/test/testRooDataSet.cxx +++ b/roofit/roofitcore/test/testRooDataSet.cxx @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -541,3 +542,73 @@ TEST(RooDataSet, ReduceCompositeDataStoreByIndexCat) EXPECT_EQ(dsCat->getCurrentIndex(), initialIndexInDataStore); EXPECT_EQ(cat->getCurrentIndex(), initialIndexInWorkspace); } + +/// Reducing, merging and appending of unbinned and binned datasets, checked +/// against event counts that are computed while deterministically filling the +/// inputs. Replaces the former stressRooFit test based on the rf402 tutorial, +/// which compared against stored reference values. +TEST(RooDataSet, ReduceMergeAppend) +{ + RooRealVar x("x", "x", -10, 10); + RooRealVar y("y", "y", 0, 40); + RooCategory c("c", "c", {{"Plus", +1}, {"Minus", -1}}); + + RooDataSet d("d", "d", {x, y, c}); + + // Fill the dataset with deterministic values and keep track of the entry + // counts that the reduce operations below are expected to yield + int nYCut = 0; // entries with y > 5.17 + int nXBinCut = 0; // entries in x bins with positive bin center (10 bins in [-10, 10], so x >= 0) + for (int i = 0; i < 1000; i++) { + x.setVal(i / 50 - 10); + y.setVal(std::sqrt(1.0 * i)); + c.setLabel((i % 2) ? "Plus" : "Minus"); + d.add({x, y, c}); + if (y.getVal() > 5.17) + ++nYCut; + if (x.getVal() >= 0) + ++nXBinCut; + } + + using namespace RooFit; + + std::unique_ptr d1{d.reduce(SelectVars({x, c}))}; + std::unique_ptr d2{d.reduce(SelectVars(y))}; + std::unique_ptr d3{d.reduce(Cut("y>5.17"))}; + std::unique_ptr d4{d.reduce(SelectVars({x, c}), Cut("y>5.17"))}; + + EXPECT_EQ(d1->numEntries(), 1000); + EXPECT_EQ(d1->get()->size(), 2u); + EXPECT_EQ(d2->numEntries(), 1000); + EXPECT_EQ(d2->get()->size(), 1u); + EXPECT_EQ(d3->numEntries(), nYCut); + EXPECT_EQ(d3->get()->size(), 3u); + EXPECT_EQ(d4->numEntries(), nYCut); + EXPECT_EQ(d4->get()->size(), 2u); + + // merge() adds datasets column-wise + static_cast(*d1).merge(static_cast(d2.get())); + EXPECT_EQ(d1->numEntries(), 1000); + EXPECT_EQ(d1->get()->size(), 3u); + + // The merged column must contain the original y values + EXPECT_DOUBLE_EQ(d1->get(999)->getRealValue("y"), std::sqrt(999.)); + + // append() adds datasets row-wise + static_cast(*d1).append(static_cast(*d3)); + EXPECT_EQ(d1->numEntries(), 1000 + nYCut); + + // Binned clone of the unbinned dataset: all entries are inside the variable + // ranges, so no events may be lost + x.setBins(10); + y.setBins(10); + RooDataHist dh("dh", "binned version of d", {x, y}, d); + EXPECT_DOUBLE_EQ(dh.sumEntries(), 1000.); + + // Reduction of the binned dataset. The cut is evaluated on the bin + // centers, so the reduced histogram contains the full content of all bins + // with positive center. + std::unique_ptr dh2{dh.reduce(SelectVars(y), Cut("x>0"))}; + EXPECT_EQ(dh2->numEntries(), 10); + EXPECT_DOUBLE_EQ(dh2->sumEntries(), nXBinCut); +} diff --git a/roofit/roofitcore/test/testRooEfficiency.cxx b/roofit/roofitcore/test/testRooEfficiency.cxx new file mode 100644 index 0000000000000..8b89a95d643da --- /dev/null +++ b/roofit/roofitcore/test/testRooEfficiency.cxx @@ -0,0 +1,160 @@ +// Tests for RooEfficiency and RooEffProd +// Authors: Jonas Rembser, CERN 09/2026 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "gtest_wrapper.h" + +#include +#include + +/// The value of the efficiency pdf normalized over the acceptance category +/// must reproduce the efficiency function for the "accept" state and its +/// complement for the "reject" state. +TEST(RooEfficiency, ValueVsEfficiencyFunction) +{ + RooRealVar x("x", "x", -10, 10); + RooRealVar a("a", "a", 0.4, 0, 1); + RooRealVar b("b", "b", 5); + RooRealVar c("c", "c", -1, -10, 10); + RooFormulaVar effFunc("effFunc", "(1-a)+a*cos((x-c)/b)", RooArgList(a, b, c, x)); + + RooCategory cut("cut", "cut", {{"accept", 1}, {"reject", 0}}); + RooEfficiency effPdf("effPdf", "effPdf", effFunc, cut, "accept"); + + RooArgSet nset{cut}; + for (double xVal : {-9.5, -5., -1., 0., 2.5, 7., 9.5}) { + x.setVal(xVal); + cut.setLabel("accept"); + EXPECT_NEAR(effPdf.getVal(&nset), effFunc.getVal(), 1e-9) << "x = " << xVal; + cut.setLabel("reject"); + EXPECT_NEAR(effPdf.getVal(&nset), 1. - effFunc.getVal(), 1e-9) << "x = " << xVal; + } +} + +/// Unbinned maximum likelihood fit of an efficiency function to a dataset +/// D(x, cut), with the shape observable as conditional observable. The fit +/// must recover the efficiency parameters used for generation. Replaces the +/// former stressRooFit test based on the rf701 tutorial. +TEST_P(RooFitEvalBackendTest, ConditionalFit1D) +{ + using namespace RooFit; + + RooRealVar x("x", "x", -10, 10); + + RooRealVar a("a", "a", 0.4, 0, 1); + RooRealVar b("b", "b", 5); + RooRealVar c("c", "c", -1, -10, 10); + RooFormulaVar effFunc("effFunc", "(1-a)+a*cos((x-c)/b)", RooArgList(a, b, c, x)); + + RooCategory cut("cut", "cut", {{"accept", 1}, {"reject", 0}}); + RooEfficiency effPdf("effPdf", "effPdf", effFunc, cut, "accept"); + + // Toy model eff(cut|x) * shape(x), only needed for generation + RooPolynomial shapePdf("shapePdf", "shapePdf", x, RooConst(-0.095)); + RooProdPdf model("model", "model", shapePdf, Conditional(effPdf, cut)); + + std::unique_ptr data{model.generate({x, cut}, 10000)}; + + std::unique_ptr res{ + effPdf.fitTo(*data, ConditionalObservables(x), Save(), PrintLevel(-1), _evalBackend)}; + + EXPECT_EQ(res->status(), 0); + expectParamNear(*res, "a", 0.4); + expectParamNear(*res, "c", -1.); +} + +/// Two-dimensional version of the conditional efficiency fit. Replaces the +/// former stressRooFit test based on the rf702 tutorial. +TEST_P(RooFitEvalBackendTest, ConditionalFit2D) +{ + using namespace RooFit; + + RooRealVar x("x", "x", -10, 10); + RooRealVar y("y", "y", -10, 10); + + RooRealVar ax("ax", "ax", 0.6, 0, 1); + RooRealVar bx("bx", "bx", 5); + RooRealVar cx("cx", "cx", -1, -10, 10); + RooRealVar ay("ay", "ay", 0.2, 0, 1); + RooRealVar by("by", "by", 5); + RooRealVar cy("cy", "cy", -1, -10, 10); + + RooFormulaVar effFunc("effFunc", "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))", + RooArgList(ax, bx, cx, x, ay, by, cy, y)); + + RooCategory cut("cut", "cut", {{"accept", 1}, {"reject", 0}}); + RooEfficiency effPdf("effPdf", "effPdf", effFunc, cut, "accept"); + + RooPolynomial shapePdfX("shapePdfX", "shapePdfX", x, RooConst(-0.095)); + RooPolynomial shapePdfY("shapePdfY", "shapePdfY", y, RooConst(+0.095)); + RooProdPdf shapePdf("shapePdf", "shapePdf", RooArgSet(shapePdfX, shapePdfY)); + RooProdPdf model("model", "model", shapePdf, Conditional(effPdf, cut)); + + std::unique_ptr data{model.generate({x, y, cut}, 10000)}; + + std::unique_ptr res{ + effPdf.fitTo(*data, ConditionalObservables(RooArgSet(x, y)), Save(), PrintLevel(-1), _evalBackend)}; + + EXPECT_EQ(res->status(), 0); + expectParamNear(*res, "ax", 0.6); + expectParamNear(*res, "cx", -1.); + expectParamNear(*res, "ay", 0.2); + expectParamNear(*res, "cy", -1.); +} + +/// A pdf multiplied with an acceptance efficiency via RooEffProd. The +/// normalized value is checked against the analytically known shape, and a fit +/// to generated data must recover the decay constant. Replaces the former +/// stressRooFit test based on the rf703 tutorial. +TEST_P(RooFitEvalBackendTest, EffProdFit) +{ + using namespace RooFit; + + if (_evalBackend.value() == EvalBackend::Value::Codegen || _evalBackend.value() == EvalBackend::Value::CodegenNoGrad) { + GTEST_SKIP() << "RooEffProd is not supported by the codegen backend yet"; + } + + RooRealVar t("t", "t", 0, 5); + RooRealVar tau("tau", "tau", -1.54, -4, -0.1); + RooExponential model("model", "model", t, tau); + + // Error function to simulate a turn-on slope + RooFormulaVar eff("eff", "0.5*(std::erf((t-1)/0.5)+1)", t); + RooEffProd modelEff("modelEff", "model with efficiency", model, eff); + + // The ratio of normalized pdf values is independent of the normalization + // integral and must match exp(tau * t) * eff(t) ratios + auto rawVal = [&](double tVal) { return std::exp(tau.getVal() * tVal) * 0.5 * (std::erf((tVal - 1.) / 0.5) + 1.); }; + RooArgSet nset{t}; + auto normVal = [&](double tVal) { + t.setVal(tVal); + return modelEff.getVal(&nset); + }; + EXPECT_NEAR(normVal(1.5) / normVal(0.8), rawVal(1.5) / rawVal(0.8), 1e-6); + EXPECT_NEAR(normVal(3.0) / normVal(1.5), rawVal(3.0) / rawVal(1.5), 1e-6); + + std::unique_ptr data{modelEff.generate(t, 10000)}; + + std::unique_ptr res{modelEff.fitTo(*data, Save(), PrintLevel(-1), _evalBackend)}; + + EXPECT_EQ(res->status(), 0); + expectParamNear(*res, "tau", -1.54); +} + +INSTANTIATE_TEST_SUITE_P(RooEfficiency, RooFitEvalBackendTest, testing::Values(ROOFIT_EVAL_BACKENDS_WITH_CODEGEN), + [](testing::TestParamInfo const ¶mInfo) { + std::stringstream ss; + ss << "EvalBackend" << std::get<0>(paramInfo.param).name(); + return ss.str(); + }); diff --git a/roofit/roofitcore/test/testRooHist.cxx b/roofit/roofitcore/test/testRooHist.cxx index c7c3908df1624..36b2d7d8e01e4 100644 --- a/roofit/roofitcore/test/testRooHist.cxx +++ b/roofit/roofitcore/test/testRooHist.cxx @@ -1,16 +1,23 @@ // Tests for the RooHist // Authors: Jonas Rembser, CERN 12/2022 +#include #include +#include +#include #include #include #include +#include #include +#include #include #include +#include + /// Check that the values returned by `RooHist::getFitRangeNEvt(double xmin, /// double xmax)` are correct also for non-uniform binning. Covers ROOT-9649. TEST(RooHist, GetFitRangeNEvtWithSubrange) @@ -48,3 +55,68 @@ TEST(RooHist, GetFitRangeNEvtWithSubrange) EXPECT_FLOAT_EQ(rooHist.getFitRangeNEvt(), nEvents); EXPECT_FLOAT_EQ(rooHist.getFitRangeNEvt(x.getMin(), x.getMax()), nEvents); } + +/// RooPlot::chiSquare(), residHist() and pullHist(), validated via the +/// pull = residual / error identity and the expected chi2 behavior for an +/// exact and a deliberately distorted model, instead of comparison with stored +/// reference plots. Replaces the former stressRooFit test based on the rf109 +/// tutorial. +TEST(RooHist, ResidualsAndPulls) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + RooRandom::randomGenerator()->SetSeed(1337); + + RooWorkspace ws; + ws.factory("Gaussian::gauss(x[-10, 10], mean[0.0], sigma[3.0, 0.1, 10.0])"); + RooRealVar &x = *ws.var("x"); + RooRealVar &sigma = *ws.var("sigma"); + RooAbsPdf &gauss = *ws.pdf("gauss"); + + std::unique_ptr data{gauss.generate(x, 10000)}; + + // Frame with the model that the data was generated from + std::unique_ptr frame1{x.frame(Bins(40))}; + data->plotOn(frame1.get(), DataError(RooAbsData::SumW2)); + gauss.plotOn(frame1.get()); + + const double chi2Good = frame1->chiSquare(); + EXPECT_GT(chi2Good, 0.2); + EXPECT_LT(chi2Good, 2.0); + + // The reduced chi2 grows when the number of fit parameters is taken into + // account in the number of degrees of freedom + EXPECT_GT(frame1->chiSquare(2), chi2Good); + + // Frame with a slightly distorted model that must yield a larger chi2 + sigma.setVal(3.15); + std::unique_ptr frame2{x.frame(Bins(40))}; + data->plotOn(frame2.get(), DataError(RooAbsData::SumW2)); + gauss.plotOn(frame2.get()); + EXPECT_GT(frame2->chiSquare(), chi2Good); + + // The pulls must be the residuals divided by the data error of the + // corresponding residual point + std::unique_ptr hresid{frame2->residHist(nullptr, nullptr, false, false)}; + std::unique_ptr hpull{frame2->pullHist(nullptr, nullptr, false)}; + ASSERT_EQ(hresid->GetN(), 40); + ASSERT_EQ(hpull->GetN(), 40); + + int nChecked = 0; + for (int i = 0; i < hresid->GetN(); ++i) { + EXPECT_DOUBLE_EQ(hpull->GetPointX(i), hresid->GetPointX(i)); + const double resid = hresid->GetPointY(i); + const double err = resid > 0 ? hresid->GetErrorYlow(i) : hresid->GetErrorYhigh(i); + if (err == 0.) { + // Empty bins have zero sum-of-weights error, and the pull of a + // zero-error point is defined to be zero + EXPECT_DOUBLE_EQ(hpull->GetPointY(i), 0.) << "point " << i; + continue; + } + EXPECT_NEAR(hpull->GetPointY(i), resid / err, 1e-9) << "point " << i; + ++nChecked; + } + // The bulk of the Gaussian sample must have populated bins + EXPECT_GT(nChecked, 30); +} diff --git a/roofit/roofitcore/test/testRooProfileLL.cxx b/roofit/roofitcore/test/testRooProfileLL.cxx new file mode 100644 index 0000000000000..a3297763718cd --- /dev/null +++ b/roofit/roofitcore/test/testRooProfileLL.cxx @@ -0,0 +1,107 @@ +// Tests for RooProfileLL +// Authors: Jonas Rembser, CERN 09/2026 + +#include +#include +#include +#include +#include +#include + +#include "gtest_wrapper.h" + +#include +#include + +/// The profile likelihood estimator minimizes the likelihood with respect to +/// all floating parameters except the parameter of interest. This is checked +/// via its defining properties: it vanishes at the best-fit point, it is +/// non-negative, it is bounded from above by the plain likelihood scan, and it +/// agrees with an explicit conditional minimization. All tolerances are at the +/// scale of the Minuit convergence criterion, because RooProfileLL subtracts +/// an absolute minimum from its own internal minimization, which can deviate +/// from the externally found minimum at that level. Replaces the former +/// stressRooFit test based on the rf605 tutorial, which compared likelihood +/// scan curves against stored references. +TEST_P(RooFitEvalBackendTest, ProfileLLMinimization) +{ + using namespace RooFit; + + // Tolerance for comparing results of independent minimizations of the same + // likelihood, at the scale of the Minuit convergence (EDM) criterion + const double minimTol = 1e-3; + + RooRealVar x("x", "x", -20, 20); + + // Model with intentionally strong correlations + RooRealVar mean("mean", "mean of g1 and g2", 0, -10, 10); + RooGaussian g1("g1", "g1", x, mean, 3.0); + RooRealVar sigmaG2("sigma_g2", "width of g2", 4, 3.0, 6.0); + RooGaussian g2("g2", "g2", x, mean, sigmaG2); + RooRealVar frac("frac", "frac", 0.5, 0.0, 1.0); + RooAddPdf model("model", "model", {g1, g2}, frac); + + std::unique_ptr data{model.generate(x, 1000)}; + + std::unique_ptr nll{model.createNLL(*data, _evalBackend)}; + + // Find the global minimum + { + RooMinimizer m(*nll); + m.setPrintLevel(-1); + m.migrad(); + } + const double minNll = nll->getVal(); + const double fracBest = frac.getVal(); + + RooArgSet params{mean, sigmaG2, frac}; + RooArgSet snapshot; + params.snapshot(snapshot); + + // Plain likelihood scan in frac, with the other parameters fixed to their + // values at the global minimum. The scan points are chosen far away from + // the best-fit value of frac (around 0.65), so that the scanned likelihood + // differences are much larger than the minimization tolerance. + const std::vector scanPoints{0.2, 0.35, 0.5, 0.8}; + std::vector plainDeltaNll; + for (double v : scanPoints) { + frac.setVal(v); + plainDeltaNll.push_back(nll->getVal() - minNll); + } + params.assign(snapshot); + + RooProfileLL pll("pll_frac", "pll_frac", *nll, frac); + + // At the best-fit point the profile likelihood is zero + frac.setVal(fracBest); + EXPECT_NEAR(pll.getVal(), 0., minimTol); + + // Away from the minimum, 0 <= profile likelihood <= plain likelihood scan + for (std::size_t i = 0; i < scanPoints.size(); ++i) { + frac.setVal(scanPoints[i]); + const double p = pll.getVal(); + EXPECT_GE(p, -minimTol) << "frac = " << scanPoints[i]; + EXPECT_LE(p, plainDeltaNll[i] + minimTol) << "frac = " << scanPoints[i]; + } + + // Cross-check one profile value against an explicit minimization with the + // parameter of interest fixed + const double fracFixed = 0.3; + frac.setVal(fracFixed); + const double profileVal = pll.getVal(); + + params.assign(snapshot); + frac.setVal(fracFixed); + frac.setConstant(true); + RooMinimizer m2(*nll); + m2.setPrintLevel(-1); + m2.migrad(); + EXPECT_NEAR(profileVal, nll->getVal() - minNll, 2. * minimTol); +} + +INSTANTIATE_TEST_SUITE_P(RooProfileLL, RooFitEvalBackendTest, testing::Values(ROOFIT_EVAL_BACKENDS), + [](testing::TestParamInfo const ¶mInfo) { + std::stringstream ss; + ss << "EvalBackend" << std::get<0>(paramInfo.param).name(); + return ss.str(); + }); diff --git a/roofit/roofitcore/test/testRooRealSumPdf.cxx b/roofit/roofitcore/test/testRooRealSumPdf.cxx new file mode 100644 index 0000000000000..d0cc0a356fe69 --- /dev/null +++ b/roofit/roofitcore/test/testRooRealSumPdf.cxx @@ -0,0 +1,94 @@ +// Tests for RooRealSumPdf +// Authors: Jonas Rembser, CERN 09/2026 + +#include +#include +#include +#include +#include +#include +#include + +#include "gtest_wrapper.h" + +#include +#include + +/// A pdf defined by a sum of real-valued amplitude components. The amplitude +/// values are checked against their analytic expressions, the normalized pdf +/// is checked via normalization-independent value ratios, and a fit to +/// generated data must recover the amplitude strength. Replaces the former +/// stressRooFit test based on the rf704 tutorial, which compared against +/// stored reference plots. +TEST_P(RooFitEvalBackendTest, RealSumPdfAmplitudeSum) +{ + RooRealVar t("t", "time", -1., 15.); + RooRealVar cosa("cosa", "cos(alpha)", -1., 1.); + + RooRealVar tau("tau", "#tau", 1.5); + RooRealVar deltaGamma("deltaGamma", "deltaGamma", 0.3); + RooFormulaVar coshG("coshGBasis", "exp(-@0/ @1)*cosh(@0*@2/2)", {t, tau, deltaGamma}); + RooFormulaVar sinhG("sinhGBasis", "exp(-@0/ @1)*sinh(@0*@2/2)", {t, tau, deltaGamma}); + + RooPolyVar poly1("poly1", "poly1", cosa, RooArgList(0.5, 0.2, 0.2), 0); + RooPolyVar poly2("poly2", "poly2", cosa, RooArgList(1.0, -0.2, 3.0), 0); + + RooProduct ampl1("ampl1", "amplitude 1", {poly1, coshG}); + RooProduct ampl2("ampl2", "amplitude 2", {poly2, sinhG}); + + RooRealVar f1("f1", "f1", 1, 0, 2); + RooRealVar f2("f2", "f2", 0.5, 0, 2); + + RooRealSumPdf pdf("pdf", "pdf", RooArgList(ampl1, ampl2), RooArgList(f1, f2)); + + // Analytic expressions for the amplitude components + auto ampl1Val = [&](double tVal, double cVal) { + return (0.5 + 0.2 * cVal + 0.2 * cVal * cVal) * std::exp(-tVal / 1.5) * std::cosh(tVal * 0.3 / 2); + }; + auto ampl2Val = [&](double tVal, double cVal) { + return (1.0 - 0.2 * cVal + 3.0 * cVal * cVal) * std::exp(-tVal / 1.5) * std::sinh(tVal * 0.3 / 2); + }; + + for (double tVal : {-0.5, 0.3, 2., 8.}) { + for (double cVal : {-0.9, 0., 0.4}) { + t.setVal(tVal); + cosa.setVal(cVal); + EXPECT_NEAR(ampl1.getVal(), ampl1Val(tVal, cVal), 1e-9) << "(t, cosa) = (" << tVal << ", " << cVal << ")"; + EXPECT_NEAR(ampl2.getVal(), ampl2Val(tVal, cVal), 1e-9) << "(t, cosa) = (" << tVal << ", " << cVal << ")"; + } + } + + // Ratios of normalized pdf values are independent of the normalization + // integral and must match the raw amplitude combination + RooArgSet nset{t, cosa}; + auto normVal = [&](double tVal, double cVal) { + t.setVal(tVal); + cosa.setVal(cVal); + return pdf.getVal(&nset); + }; + auto rawVal = [&](double tVal, double cVal) { + return f1.getVal() * ampl1Val(tVal, cVal) + f2.getVal() * ampl2Val(tVal, cVal); + }; + EXPECT_NEAR(normVal(2., 0.4) / normVal(1., -0.5), rawVal(2., 0.4) / rawVal(1., -0.5), 1e-6); + EXPECT_NEAR(normVal(6., 0.8) / normVal(2., 0.4), rawVal(6., 0.8) / rawVal(2., 0.4), 1e-6); + + // Generate toy data and fit with one amplitude strength floating. Only the + // relative amplitude strength is defined, so f1 is kept constant. + std::unique_ptr data{pdf.generate({t, cosa}, 10000)}; + + f1.setConstant(true); + f2.setVal(1.0); + + std::unique_ptr res{pdf.fitTo(*data, RooFit::Save(), RooFit::PrintLevel(-1), _evalBackend)}; + + EXPECT_EQ(res->status(), 0); + EXPECT_EQ(res->covQual(), 3); + expectParamNear(*res, "f2", 0.5); +} + +INSTANTIATE_TEST_SUITE_P(RooRealSumPdf, RooFitEvalBackendTest, testing::Values(ROOFIT_EVAL_BACKENDS), + [](testing::TestParamInfo const ¶mInfo) { + std::stringstream ss; + ss << "EvalBackend" << std::get<0>(paramInfo.param).name(); + return ss.str(); + });