[RF][RS] Store fit results in the AsymptoticCalculator

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

Author: guitargeek

roofit/roostats/inc/RooStats/AsymptoticCalculator.h

@@ -15,8 +15,11 @@
 #include "RooArgSet.h"
 #include "Rtypes.h"
 
+#include <memory>
+
 class RooArgList;
 class RooCategory;
+class RooFitResult;
 class RooRealVar;
 class RooPoisson;
 class RooProdPdf;
@@ -100,6 +103,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:
@@ -117,6 +143,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<RooFitResult> fFitResultUncondObs;    ///<! result of unconditional fit to observed data
+      mutable std::unique_ptr<RooFitResult> fFitResultCondObs;      ///<! result of conditional fit to observed data
+      mutable std::unique_ptr<RooFitResult> fFitResultUncondAsimov; ///<! result of fit to Asimov data at the alt POI
+      mutable std::unique_ptr<RooFitResult> fFitResultCondAsimov;   ///<! result of conditional fit to Asimov data
+
       ClassDefOverride(AsymptoticCalculator,0)
    };
 }

roofit/roostats/src/AsymptoticCalculator.cxx

@@ -35,8 +35,18 @@ If more than one POI exists, only the first one is used.
 The calculator can generate Asimov datasets from two kinds of PDFs:
 - "Counting" distributions: RooPoisson, RooGaussian, or products of RooPoissons.
 - Extended, *i.e.* number of events can be read off from extended likelihood term.
-*/
 
+The fits performed by the calculator can be steered with the global default
+minimizer options, *e.g.* via ROOT::Math::MinimizerOptions::SetDefaultStrategy()
+and ROOT::Math::MinimizerOptions::SetDefaultTolerance() (the tolerance is
+clamped to a minimum value of 1). The RooFitResult objects of the fits are
+retrievable after calling GetHypoTest() via GetFitResultUncondObs(),
+GetFitResultCondObs(), GetFitResultUncondAsimov() and GetFitResultCondAsimov(),
+so quantities like the minimizer status or the EDM at the minimum can be
+inspected, for example to cross-check a fit that did not converge. Note that
+when the calculator is driven by the HypoTestInverter, the stored conditional
+fit results correspond to the last scanned point.
+*/
 
 #include "RooStats/AsymptoticCalculator.h"
 #include "RooStats/ModelConfig.h"
@@ -84,7 +94,8 @@ int &fgPrintLevel()
 }
 
 // Forward declaration.
-double EvaluateNLL(RooStats::ModelConfig const &modelConfig, RooAbsData &data, const RooArgSet *poiSet = nullptr);
+double EvaluateNLL(RooStats::ModelConfig const &modelConfig, RooAbsData &data, const RooArgSet *poiSet = nullptr,
+                   std::unique_ptr<RooFitResult> *fitResult = nullptr);
 
 } // namespace
 
@@ -196,11 +207,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<RooRealVar*>(fBestFitPoi.first());
@@ -283,7 +302,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");
@@ -300,10 +323,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<RooFitResult> *fitResult)
 {
     int verbose = fgPrintLevel();
 
+    if (fitResult)
+       fitResult->reset();
+
     RooAbsPdf &pdf = *modelConfig.GetPdf();
 
     RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
@@ -415,13 +442,12 @@ double EvaluateNLL(RooStats::ModelConfig const& modelConfig, RooAbsData& data, c
           }
        }
 
-       std::unique_ptr<RooFitResult> 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<RooFitResult> 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<RooFitResult>{minim.save()};
-       }
-       if (result){
+       if (status >= 0 && result) {
           if (RooStats::NLLOffsetMode() != "initial") {
              val = result->minNll();
           } else {
@@ -431,12 +457,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);
     }
 
@@ -528,7 +556,12 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {
    }
 
    // evaluate the conditional NLL on the observed data for the snapshot value
-   double condNLL = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData&>(*GetData()), &poiTest);
+   double condNLL = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData &>(*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);
 
@@ -550,14 +583,20 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {
             << "AsymptoticCalculator:  unconditional fit failed before - retry to do it now " << std::endl;
       }
 
-      double nll = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData&>(*GetData()));
+      std::unique_ptr<RooFitResult> refitResult;
+      double nll = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData &>(*GetData()), nullptr, &refitResult);
 
       if (nll < fNLLObs || (TMath::IsNaN(fNLLObs) && !TMath::IsNaN(nll) ) ) {
          oocoutW(nullptr,Minimization) << "AsymptoticCalculator:  Found a better unconditional minimum "
                                            << " old NLL = " << fNLLObs << " old muHat " << muHat->getVal() << std::endl;
 
          // 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();
@@ -612,8 +651,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  );
 
@@ -632,14 +675,20 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {
            << std::endl;
       }
 
-      double nll = EvaluateNLL(*GetNullModel(), *fAsimovData);
+      std::unique_ptr<RooFitResult> 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"
                                            << " old NLL = " << fNLLAsimov << std::endl;
 
          // 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;

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

roofit/roostats/test/testAsymptoticCalculator.cxx

@@ -1,6 +1,7 @@
 // Author: Jonas Rembser, CERN  01/2025
 
 #include "RooDataSet.h"
+#include "RooFitResult.h"
 #include "RooMultiVarGaussian.h"
 #include "RooRealVar.h"
 #include "RooWorkspace.h"
@@ -9,6 +10,7 @@
 #include "RooStats/ModelConfig.h"
 
 #include "Math/ProbFuncMathCore.h"
+#include "ROOT/TestSupport.hxx"
 
 #include "gtest/gtest.h"
 
@@ -147,3 +149,83 @@ TEST(AsymptoticCalculator, CountingAsimovDataSetFloatingParams)
    ws.var("mean")->setConstant(true);
    checkAsimov("gauss1", 20.0);
 }
+
+// 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<ModelConfig> bModel{static_cast<ModelConfig *>(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<HypoTestResult> 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<RooRealVar *>(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<RooRealVar *>(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);
+}