diff --git a/roofit/roofitcore/inc/RooExtendPdf.h b/roofit/roofitcore/inc/RooExtendPdf.h index 21706c59dd46c..09bcd1f2f3b71 100644 --- a/roofit/roofitcore/inc/RooExtendPdf.h +++ b/roofit/roofitcore/inc/RooExtendPdf.h @@ -48,6 +48,8 @@ class RooExtendPdf : public RooAbsPdf { double expectedEvents(const RooArgSet* nset) const override ; std::unique_ptr createExpectedEventsFunc(const RooArgSet* nset) const override; + std::unique_ptr compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override; + RooAbsPdf const& pdf() const { return *_pdf; } RooAbsReal const &getN() const { return *_n; } diff --git a/roofit/roofitcore/src/RooExtendPdf.cxx b/roofit/roofitcore/src/RooExtendPdf.cxx index b45f7e3102016..60ebfada4c233 100644 --- a/roofit/roofitcore/src/RooExtendPdf.cxx +++ b/roofit/roofitcore/src/RooExtendPdf.cxx @@ -49,6 +49,8 @@ the nominal integration range \f$ \mathrm{normRegion}[x] \f$. #include "RooFitImplHelpers.h" +#include + using std::endl; @@ -136,6 +138,47 @@ double RooExtendPdf::expectedEvents(const RooArgSet* nset) const } +/// RooExtendPdf is self-normalized (selfNormalized() returns true), which means +/// it delegates the normalization of the shape entirely to the wrapped pdf. +/// The fractional correction of the extended term in expectedEvents() / +/// createExpectedEventsFunc() is also evaluated from the wrapped pdf. For a +/// likelihood over a sub-range this means the wrapped pdf must be normalized +/// over the same range as this RooExtendPdf (the fit range), otherwise both the +/// shape term (biasing the shape parameters) and the extended term (collapsing +/// the yield to the fit-range count instead of reinterpreting it to the range +/// requested in the constructor) are wrong. +/// +/// Since the normalization range of a pdf is not automatically propagated to +/// its servers, we propagate it explicitly to the wrapped pdf while the +/// computation graph is compiled, mirroring what RooAddPdf does for its +/// component pdfs. This is only done in likelihood mode: for a binned chi2, the +/// per-bin predictions must stay normalized over the full range so that the +/// chi2 remains additive across sub-ranges. See GitHub issue #22959. +std::unique_ptr +RooExtendPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext &ctx) const +{ + // RAII to temporarily set the wrapped pdf's normalization range to the one + // of this RooExtendPdf and reset it after compilation. + struct NormRangeGuard { + NormRangeGuard(RooAbsPdf &pdf, const char *rangeName) + : _pdf{pdf}, _oldNormRange{pdf.normRange() ? pdf.normRange() : ""} + { + _pdf.setNormRange(rangeName); + } + ~NormRangeGuard() { _pdf.setNormRange(_oldNormRange.empty() ? nullptr : _oldNormRange.c_str()); } + RooAbsPdf &_pdf; + std::string _oldNormRange; + }; + + if (!ctx.likelihoodMode()) { + return RooAbsPdf::compileForNormSet(normSet, ctx); + } + + NormRangeGuard guard{*_pdf, normRange()}; + return RooAbsPdf::compileForNormSet(normSet, ctx); +} + + std::unique_ptr RooExtendPdf::createExpectedEventsFunc(const RooArgSet *nset) const { const RooAbsPdf& pdf = *_pdf; diff --git a/roofit/roofitcore/test/testTestStatistics.cxx b/roofit/roofitcore/test/testTestStatistics.cxx index d663fd6601886..ef47f76d0a501 100644 --- a/roofit/roofitcore/test/testTestStatistics.cxx +++ b/roofit/roofitcore/test/testTestStatistics.cxx @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -1212,3 +1213,61 @@ TEST(RooChi2Var, BinnedRangeAdditivityAndNormalization) EXPECT_NEAR(chi2ExtHi, chi2FuncHi, tol); EXPECT_NEAR(chi2ExtLoHi, chi2FuncLoHi, tol); } + +// Extended likelihood fit of a RooExtendPdf in a sub-range. This is the +// scenario of the rf204b_extendedLikelihood_rangedFit tutorial. When the +// RooExtendPdf is constructed with a range name, the fitted yield must be +// reinterpreted to that range (here the full range) instead of collapsing to +// the number of events observed in the fit range, and the shape parameters must +// not be biased. Both must agree with the equivalent RooAddPdf model, whose +// coefficients are automatically reinterpreted to the full range. Regression +// test for GitHub issue #22959. +TEST(RooNLLVar, RooExtendPdfRangedFit) +{ + using namespace RooFit; + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRandom::randomGenerator()->SetSeed(42); + + RooRealVar x("x", "x", 10, 100); + x.setRange("LEFT", 10, 20); + x.setRange("RIGHT", 60, 100); + x.setRange("FULL", 10, 100); + + RooRealVar alphaGen("alpha", "alpha", -0.04, -0.1, -0.0); + RooGenericPdf modelGen("model", "exp(alpha*x)", {x, alphaGen}); + std::unique_ptr data{modelGen.generate(x, 10000)}; + + const double nFull = data->sumEntries(); + const double nFitRange = std::unique_ptr{data->reduce(CutRange("LEFT,RIGHT"))}->sumEntries(); + // The fit range only contains a fraction of the events, so this test is only + // meaningful if the two counts are clearly different. + ASSERT_LT(nFitRange, 0.8 * nFull); + + // Extended fit with RooExtendPdf, yield N interpreted over the FULL range. + RooRealVar alpha1("alpha", "alpha", -0.04, -0.1, -0.0); + RooGenericPdf model1("model", "exp(alpha*x)", {x, alpha1}); + RooRealVar nExpected("N", "Extended term", 5000, 0, 20000); + RooExtendPdf extmodel("extmodel", "Extended model", model1, nExpected, "FULL"); + std::unique_ptr resExt{extmodel.fitTo(*data, Range("LEFT,RIGHT"), PrintLevel(-1), Save())}; + + // Equivalent extended fit with a single-component RooAddPdf, whose + // coefficient is automatically reinterpreted to the full range in a ranged + // fit. + RooRealVar alpha2("alpha", "alpha", -0.04, -0.1, -0.0); + RooGenericPdf model2("model", "exp(alpha*x)", {x, alpha2}); + RooRealVar nBkg("Nbkg", "Number of background events", 5000, 0, 20000); + RooAddPdf modelsum("modelsum", "", RooArgList{model2}, RooArgList{nBkg}); + std::unique_ptr resAdd{modelsum.fitTo(*data, Range("LEFT,RIGHT"), PrintLevel(-1), Save())}; + + // The RooExtendPdf yield must match the RooAddPdf background yield... + EXPECT_NEAR(nExpected.getVal(), nBkg.getVal(), 1e-3 * nBkg.getVal()); + // ...and both refer to the full range, not the (much smaller) fit-range count. + EXPECT_NEAR(nExpected.getVal(), nFull, 0.05 * nFull); + EXPECT_GT(nExpected.getVal(), 1.5 * nFitRange); + + // The shape parameter must be unbiased, i.e. agree with the RooAddPdf fit and + // with the generated value. + EXPECT_NEAR(alpha1.getVal(), alpha2.getVal(), 1e-4); + EXPECT_NEAR(alpha1.getVal(), -0.04, 5. * alpha1.getError()); +} diff --git a/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.C b/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.C index c826fe1b2b83a..9182b82ba7381 100644 --- a/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.C +++ b/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.C @@ -83,8 +83,8 @@ void rf204b_extendedLikelihood_rangedFit() // If the fitted domain is a subset of `FULL`, though, the expected event count is divided by // \f[ // \mathrm{frac} = \frac{ - // \int_{\mathrm{Fit range}} \mathrm{model}(x) \; \mathrm{d}x }{ - // \int_{\mathrm{Full range}} \mathrm{model}(x) \; \mathrm{d}x }. + // \int_{\mathrm{Full range}} \mathrm{model}(x) \; \mathrm{d}x }{ + // \int_{\mathrm{Fit range}} \mathrm{model}(x) \; \mathrm{d}x }. // \f] // `N` will therefore return the count extrapolated to the full range instead of the fit range. // @@ -159,10 +159,9 @@ void rf204b_extendedLikelihood_rangedFit() // This model will automatically insert the correction factor for the reinterpretation of Nsig and Nnbkg in the full ranges. // - // When this happens, it reports this with lines like the following: + // When this happens, it reports this with a line like the following: // ``` - // [#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_LEFT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables - // [#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_RIGHT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables + // [#1] INFO:Fitting -- RooAbsPdf::fitTo(modelsum) fixing normalization set for coefficient determination to observables in data // ``` std::unique_ptr r3{modelsum.fitTo(*data, Range("LEFT,RIGHT"), PrintLevel(-1), Save())}; diff --git a/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.py b/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.py index 07d2c80b4fcbd..19b6d55a12dcc 100644 --- a/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.py +++ b/tutorials/roofit/roofit/rf204b_extendedLikelihood_rangedFit.py @@ -72,8 +72,8 @@ # If the fitted domain is a subset of `FULL`, though, the expected event count is divided by # \f[ # \mathrm{frac} = \frac{ -# \int_{\mathrm{Fit range}} \mathrm{model}(x) \; \mathrm{d}x }{ -# \int_{\mathrm{Full range}} \mathrm{model}(x) \; \mathrm{d}x }. +# \int_{\mathrm{Full range}} \mathrm{model}(x) \; \mathrm{d}x }{ +# \int_{\mathrm{Fit range}} \mathrm{model}(x) \; \mathrm{d}x }. # \f] # `N` will therefore return the count extrapolated to the full range instead of the fit range. # @@ -152,10 +152,9 @@ # This model will automatically insert the correction factor for the # reinterpretation of Nsig and Nnbkg in the full ranges. # -# When this happens, it reports this with lines like the following: +# When this happens, it reports this with a line like the following: # ``` -# [#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_LEFT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables -# [#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_RIGHT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables +# [#1] INFO:Fitting -- RooAbsPdf::fitTo(modelsum) fixing normalization set for coefficient determination to observables in data # ``` r3 = modelsum.fitTo(data, Range="LEFT,RIGHT", PrintLevel=-1, Save=True)