From 66a3875442caf40163a5f05e0310e8218e3f82f7 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 19 Aug 2026 11:37:30 +0000 Subject: [PATCH] [RF][RS] Store fit results in the AsymptoticCalculator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The AsymptoticCalculator performed all its fits in a local helper that saved a RooFitResult only to read off the minimum NLL value and then discarded it. Users therefore had no way to check programmatically whether the underlying fits actually converged, even though the calculator happily reports a significance also when they did not. Keep the RooFitResult objects of the four fits (unconditional and conditional, on observed and Asimov data) as members of the calculator and expose them with the new getters GetFitResultUncondObs(), GetFitResultCondObs(), GetFitResultUncondAsimov() and GetFitResultCondAsimov(). The results are now saved also when the minimization failed, so that the minimizer status and EDM of a non-converged fit can be inspected. The stored unconditional results are updated when GetHypoTest() finds a better minimum in its refit fallback. Also document in the class description that the fits can already be steered via the ROOT::Math::MinimizerOptions defaults for strategy and tolerance, which was the second request in the ticket. Closes JIRA [ROOT-10066](https://its.cern.ch/jira/browse/ROOT-10066). 🤖 Done with the help of AI --- .../inc/RooStats/AsymptoticCalculator.h | 31 +++++++ roofit/roostats/src/AsymptoticCalculator.cxx | 83 +++++++++++++---- roofit/roostats/test/CMakeLists.txt | 2 +- .../test/testAsymptoticCalculator.cxx | 88 +++++++++++++++++++ 4 files changed, 186 insertions(+), 18 deletions(-) diff --git a/roofit/roostats/inc/RooStats/AsymptoticCalculator.h b/roofit/roostats/inc/RooStats/AsymptoticCalculator.h index fd0b5e0477479..fb9e20b7ba06a 100644 --- a/roofit/roostats/inc/RooStats/AsymptoticCalculator.h +++ b/roofit/roostats/inc/RooStats/AsymptoticCalculator.h @@ -15,8 +15,11 @@ #include "RooArgSet.h" #include "Rtypes.h" +#include + class RooArgList; class RooCategory; +class RooFitResult; class RooRealVar; class RooPoisson; class RooProdPdf; @@ -93,6 +96,29 @@ namespace RooStats { /// return best fit value for all parameters const RooArgSet & GetBestFitParams() const { return fBestFitPoi; } + /// Result of the unconditional fit to the observed data, performed by + /// Initialize() (updated if GetHypoTest() finds a better minimum). + /// Returns nullptr if the fit was skipped or has not been run yet. + /// The calculator keeps ownership of the returned object. + const RooFitResult *GetFitResultUncondObs() const { return fFitResultUncondObs.get(); } + /// Result of the conditional fit to the observed data with the POI fixed + /// to the tested value, from the last call to GetHypoTest(). + /// Returns nullptr if the fit was skipped or has not been run yet. + /// The calculator keeps ownership of the returned object. + const RooFitResult *GetFitResultCondObs() const { return fFitResultCondObs.get(); } + /// Result of the fit to the Asimov data set with the POI fixed to the + /// value of the alternate-model snapshot, performed by Initialize(). Since + /// the Asimov data set is generated at that POI value, this corresponds to + /// the unconditional minimum (updated if GetHypoTest() finds a better + /// minimum). Returns nullptr if the fit was skipped or has not been run + /// yet. The calculator keeps ownership of the returned object. + const RooFitResult *GetFitResultUncondAsimov() const { return fFitResultUncondAsimov.get(); } + /// Result of the conditional fit to the Asimov data set with the POI + /// fixed to the tested value, from the last call to GetHypoTest(). + /// Returns nullptr if the fit was skipped or has not been run yet. + /// The calculator keeps ownership of the returned object. + const RooFitResult *GetFitResultCondAsimov() const { return fFitResultCondAsimov.get(); } + static void SetPrintLevel(int level); private: @@ -109,6 +135,11 @@ namespace RooStats { mutable RooArgSet fBestFitPoi; ///< snapshot of best fitted POI values mutable RooArgSet fBestFitParams; ///< snapshot of all best fitted Parameter values + mutable std::unique_ptr fFitResultUncondObs; /// fFitResultCondObs; /// fFitResultUncondAsimov; /// fFitResultCondAsimov; /// *fitResult = nullptr); } // namespace @@ -198,11 +209,19 @@ bool AsymptoticCalculator::Initialize() const { fBestFitPoi.removeAll(); fBestFitParams.removeAll(); fAsimovGlobObs.removeAll(); + fFitResultUncondObs.reset(); + fFitResultCondObs.reset(); + fFitResultUncondAsimov.reset(); + fFitResultCondAsimov.reset(); // evaluate the unconditional nll for the full model on the observed data if (verbose >= 0) oocoutP(nullptr,Eval) << "AsymptoticCalculator::Initialize - Find best unconditional NLL on observed data" << std::endl; - fNLLObs = EvaluateNLL(*GetNullModel(), data); + fNLLObs = EvaluateNLL(*GetNullModel(), data, nullptr, &fFitResultUncondObs); + if (fFitResultUncondObs) { + fFitResultUncondObs->SetName("fitResultUncondObs"); + fFitResultUncondObs->SetTitle("Unconditional fit to observed data"); + } // fill also snapshot of best poi poi->snapshot(fBestFitPoi); RooRealVar * muBest = dynamic_cast(fBestFitPoi.first()); @@ -285,7 +304,11 @@ bool AsymptoticCalculator::Initialize() const { << muAlt->GetName() << " ) = " << muAlt->getVal() << std::endl; } - fNLLAsimov = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiAlt ); + fNLLAsimov = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiAlt, &fFitResultUncondAsimov); + if (fFitResultUncondAsimov) { + fFitResultUncondAsimov->SetName("fitResultUncondAsimov"); + fFitResultUncondAsimov->SetTitle("Fit to Asimov data with POI fixed to the alt-model snapshot"); + } // for unconditional fit //fNLLAsimov = EvaluateNLL( *nullPdf, *fAsimovData); //poi->Print("v"); @@ -302,10 +325,14 @@ bool AsymptoticCalculator::Initialize() const { namespace { -double EvaluateNLL(RooStats::ModelConfig const& modelConfig, RooAbsData& data, const RooArgSet *poiSet) +double EvaluateNLL(RooStats::ModelConfig const &modelConfig, RooAbsData &data, const RooArgSet *poiSet, + std::unique_ptr *fitResult) { int verbose = fgPrintLevel(); + if (fitResult) + fitResult->reset(); + RooAbsPdf &pdf = *modelConfig.GetPdf(); RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow(); @@ -417,13 +444,12 @@ double EvaluateNLL(RooStats::ModelConfig const& modelConfig, RooAbsData& data, c } } - std::unique_ptr result; + // save the fit result also in case of failure, so that the status of a + // non-converged fit can be inspected by the user + std::unique_ptr result{minim.save()}; // ignore errors in Hesse or in Improve and also when matrix was made pos def (status returned = 1) - if (status >= 0) { - result = std::unique_ptr{minim.save()}; - } - if (result){ + if (status >= 0 && result) { if (RooStats::NLLOffsetMode() != "initial") { val = result->minNll(); } else { @@ -433,12 +459,14 @@ double EvaluateNLL(RooStats::ModelConfig const& modelConfig, RooAbsData& data, c if (!previous) RooAbsReal::setHideOffset(false) ; } - } - else { + } else { oocoutE(nullptr,Fitting) << "FIT FAILED !- return a NaN NLL " << std::endl; val = TMath::QuietNaN(); } + if (fitResult) + *fitResult = std::move(result); + minim.optimizeConst(false); } @@ -530,7 +558,12 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { } // evaluate the conditional NLL on the observed data for the snapshot value - double condNLL = EvaluateNLL(*GetNullModel(), const_cast(*GetData()), &poiTest); + double condNLL = EvaluateNLL(*GetNullModel(), const_cast(*GetData()), &poiTest, &fFitResultCondObs); + if (fFitResultCondObs) { + fFitResultCondObs->SetName("fitResultCondObs"); + fFitResultCondObs->SetTitle( + TString::Format("Conditional fit to observed data for %s = %g", muTest->GetName(), muTest->getVal())); + } double qmu = 2.*(condNLL - fNLLObs); @@ -552,7 +585,8 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { << "AsymptoticCalculator: unconditional fit failed before - retry to do it now " << std::endl; } - double nll = EvaluateNLL(*GetNullModel(), const_cast(*GetData())); + std::unique_ptr refitResult; + double nll = EvaluateNLL(*GetNullModel(), const_cast(*GetData()), nullptr, &refitResult); if (nll < fNLLObs || (TMath::IsNaN(fNLLObs) && !TMath::IsNaN(nll) ) ) { oocoutW(nullptr,Minimization) << "AsymptoticCalculator: Found a better unconditional minimum " @@ -560,6 +594,11 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { // update values fNLLObs = nll; + if (refitResult) { + fFitResultUncondObs = std::move(refitResult); + fFitResultUncondObs->SetName("fitResultUncondObs"); + fFitResultUncondObs->SetTitle("Unconditional fit to observed data"); + } const RooArgSet * poi = GetNullModel()->GetParametersOfInterest(); assert(poi); fBestFitPoi.removeAll(); @@ -614,8 +653,12 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { if (verbose > 0) oocoutP(nullptr,Eval) << "AsymptoticCalculator::GetHypoTest -- Find best conditional NLL on ASIMOV data set .... " << std::endl; - double condNLL_A = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiTest); - + double condNLL_A = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiTest, &fFitResultCondAsimov); + if (fFitResultCondAsimov) { + fFitResultCondAsimov->SetName("fitResultCondAsimov"); + fFitResultCondAsimov->SetTitle( + TString::Format("Conditional fit to Asimov data for %s = %g", muTest->GetName(), muTest->getVal())); + } double qmu_A = 2.*(condNLL_A - fNLLAsimov ); @@ -634,7 +677,8 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { << std::endl; } - double nll = EvaluateNLL(*GetNullModel(), *fAsimovData); + std::unique_ptr refitResult; + double nll = EvaluateNLL(*GetNullModel(), *fAsimovData, nullptr, &refitResult); if (nll < fNLLAsimov || (TMath::IsNaN(fNLLAsimov) && !TMath::IsNaN(nll) )) { oocoutW(nullptr,Minimization) << "AsymptoticCalculator: Found a better unconditional minimum for Asimov data set" @@ -642,6 +686,11 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { // update values fNLLAsimov = nll; + if (refitResult) { + fFitResultUncondAsimov = std::move(refitResult); + fFitResultUncondAsimov->SetName("fitResultUncondAsimov"); + fFitResultUncondAsimov->SetTitle("Unconditional fit to Asimov data"); + } oocoutW(nullptr,Minimization) << "AsymptoticCalculator: New minimum found for " << " NLL = " << fNLLAsimov << std::endl; diff --git a/roofit/roostats/test/CMakeLists.txt b/roofit/roostats/test/CMakeLists.txt index 7770dc075e93e..e14c02d8a35e6 100644 --- a/roofit/roostats/test/CMakeLists.txt +++ b/roofit/roostats/test/CMakeLists.txt @@ -1,4 +1,4 @@ -ROOT_ADD_GTEST(testAsymptoticCalculator testAsymptoticCalculator.cxx LIBRARIES RooStats) +ROOT_ADD_GTEST(testAsymptoticCalculator testAsymptoticCalculator.cxx LIBRARIES RooStats ROOT::TestSupport) ROOT_ADD_GTEST(testBayesianCalculator testBayesianCalculator.cxx LIBRARIES RooStats) ROOT_ADD_GTEST(testHypoTestInvResult testHypoTestInvResult.cxx LIBRARIES RooStats diff --git a/roofit/roostats/test/testAsymptoticCalculator.cxx b/roofit/roostats/test/testAsymptoticCalculator.cxx index 6d97b7daa5c3b..53d9aa2eb90fd 100644 --- a/roofit/roostats/test/testAsymptoticCalculator.cxx +++ b/roofit/roostats/test/testAsymptoticCalculator.cxx @@ -1,7 +1,15 @@ // Author: Jonas Rembser, CERN 01/2025 +#include "RooDataSet.h" +#include "RooFitResult.h" #include "RooMultiVarGaussian.h" +#include "RooRealVar.h" +#include "RooWorkspace.h" #include "RooStats/AsymptoticCalculator.h" +#include "RooStats/HypoTestResult.h" +#include "RooStats/ModelConfig.h" + +#include "ROOT/TestSupport.hxx" #include "gtest/gtest.h" @@ -37,3 +45,83 @@ TEST(AsymptoticCalculator, CountingAsimovDataSetFromMultiVarGaussian) EXPECT_EQ(dataX.getVal(), mu.getVal()); } } + +// Check that the fit results of the fits performed in Initialize() and +// GetHypoTest() are stored and can be accessed by the user (JIRA ROOT-10066). +TEST(AsymptoticCalculator, StoredFitResults) +{ + using namespace RooStats; + + // On/off Poisson counting model: signal region with s + b expected events, + // control region constraining the background via tau * b. + RooWorkspace ws; + ws.factory("Poisson::px(x[150, 0, 500], sum::splusb(s[0, 0, 100], b[100, 0, 300]))"); + ws.factory("Poisson::py(y[100, 0, 500], prod::taub(tau[1.0], b))"); + ws.factory("PROD::model(px, py)"); + + RooRealVar &s = *ws.var("s"); + RooArgSet obs{*ws.var("x"), *ws.var("y")}; + + RooDataSet data{"data", "data", obs}; + data.add(obs); + + ModelConfig sbModel{"sbModel", &ws}; + sbModel.SetPdf(*ws.pdf("model")); + sbModel.SetObservables(obs); + sbModel.SetParametersOfInterest(RooArgSet{s}); + sbModel.SetNuisanceParameters(RooArgSet{*ws.var("b")}); + s.setVal(50.0); + sbModel.SetSnapshot(RooArgSet{s}); + + std::unique_ptr bModel{static_cast(sbModel.Clone("bModel"))}; + s.setVal(0.0); + bModel->SetSnapshot(RooArgSet{s}); + + // Some of the fits start at parameter values that already correspond to the + // minimum, in which case Minuit2 emits a harmless line-search warning. + ROOT::TestSupport::CheckDiagsRAII checkDiag; + checkDiag.optionalDiag(kWarning, "Minuit2", "VariableMetricBuilder No improvement in line search", false); + + AsymptoticCalculator calc{data, *bModel, sbModel}; + calc.SetOneSided(true); + + // Before running the hypothesis test, no fit results are available. + EXPECT_EQ(calc.GetFitResultCondObs(), nullptr); + EXPECT_EQ(calc.GetFitResultCondAsimov(), nullptr); + + std::unique_ptr result{calc.GetHypoTest()}; + ASSERT_NE(result, nullptr); + + const RooFitResult *uncondObs = calc.GetFitResultUncondObs(); + const RooFitResult *condObs = calc.GetFitResultCondObs(); + const RooFitResult *uncondAsimov = calc.GetFitResultUncondAsimov(); + const RooFitResult *condAsimov = calc.GetFitResultCondAsimov(); + + ASSERT_NE(uncondObs, nullptr); + ASSERT_NE(condObs, nullptr); + ASSERT_NE(uncondAsimov, nullptr); + ASSERT_NE(condAsimov, nullptr); + + EXPECT_EQ(uncondObs->status(), 0); + EXPECT_EQ(condObs->status(), 0); + EXPECT_EQ(uncondAsimov->status(), 0); + EXPECT_EQ(condAsimov->status(), 0); + + // The best-fit POI value from the unconditional fit result must be + // consistent with the one stored by the calculator. + auto *sFitUncond = static_cast(uncondObs->floatParsFinal().find("s")); + ASSERT_NE(sFitUncond, nullptr); + EXPECT_DOUBLE_EQ(sFitUncond->getVal(), calc.GetMuHat()->getVal()); + + // In the conditional fits the POI is fixed to the tested value from the + // null-model snapshot, so it appears in the constant parameter list. + auto *sFitCond = static_cast(condObs->constPars().find("s")); + ASSERT_NE(sFitCond, nullptr); + EXPECT_DOUBLE_EQ(sFitCond->getVal(), 50.0); + + // The profile likelihood ratio test statistic reconstructed from the stored + // fit results must be non-negative, up to the same numerical tolerance that + // the calculator itself uses for qmu. + EXPECT_GE(condObs->minNll() - uncondObs->minNll(), -1.e-3); + EXPECT_GE(condAsimov->minNll() - uncondAsimov->minNll(), -1.e-3); +}