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(); + });