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