[RF] Fix counting Asimov generation with multiple floating parameters

When generating a counting Asimov dataset, the AsymptoticCalculator
identified the expected value of an observable in a Gaussian or Poisson
term as "the single non-constant server that is not the observable".
This heuristic failed whenever both the mean and the width of a Gaussian
were floating, could not see constness through derived quantities like
RooFormulaVars, and silently set the observable to the value of the
*sigma* parameter when the mean was constant but sigma floated.

The heuristic is unnecessary: RooGaussian and RooPoisson expose their
proxies via getX() and getMean(). Use the exact (x, mean) pair and set
whichever of the two is the observable to the value of the other one.
The width and the constness of the parameters no longer matter, and the
direction where the mean is the observable (constraint terms with global
observables) keeps working. The RooMultiVarGaussian path goes through
the same helper with xVec()[i] and muVec()[i].

Fixes the four failure cases from ROOT-10096, covered by a new test.

Fixes [ROOT-10096](https://its.cern.ch/jira/browse/ROOT-10069)

🤖 Done with the help of AI

Author: guitargeek

roofit/roostats/src/AsymptoticCalculator.cxx

@@ -69,8 +69,6 @@ The calculator can generate Asimov datasets from two kinds of PDFs:
 
 #include "TStopwatch.h"
 
-#include <ROOT/RSpan.hxx>
-
 using namespace RooStats;
 using std::string, std::unique_ptr;
 
@@ -939,41 +937,31 @@ void FillBins(const RooAbsPdf & pdf, const RooArgList &obs, RooAbsData & data, i
 
 }
 
-bool setObsToExpected(std::span<RooAbsArg *> servers, const RooArgSet &obs, std::string const &errPrefix)
+bool setObsToExpected(RooAbsArg &x, RooAbsArg &mean, const RooArgSet &obs, std::string const &errPrefix)
 {
-   RooRealVar *myobs = nullptr;
-   RooAbsReal *myexp = nullptr;
-   for (RooAbsArg *a : servers) {
-      if (obs.contains(*a)) {
-         if (myobs != nullptr) {
-            oocoutF(nullptr,Generation) << errPrefix << "Has two observables ?? " << std::endl;
-            return false;
-         }
-         myobs = dynamic_cast<RooRealVar *>(a);
-         if (myobs == nullptr) {
-            oocoutF(nullptr,Generation) << errPrefix << "Observable is not a RooRealVar??" << std::endl;
-            return false;
-         }
-      } else {
-         if (!a->isConstant() ) {
-            if (myexp != nullptr) {
-               oocoutE(nullptr,Generation) << errPrefix << "Has two non-const arguments  " << std::endl;
-               return false;
-            }
-            myexp = dynamic_cast<RooAbsReal *>(a);
-            if (myexp == nullptr) {
-               oocoutF(nullptr,Generation) << errPrefix << "Expected is not a RooAbsReal??" << std::endl;
-               return false;
-            }
-         }
-      }
+   // Figure out which of the two arguments is the observable that should be
+   // set to the expected value given by the other one. Usually the observable
+   // is "x", but also the mean parameter can be the observable: this happens
+   // for example in constraint terms, where the global observable takes the
+   // role of the mean.
+   const bool xIsObs = obs.contains(x);
+   const bool meanIsObs = obs.contains(mean);
+   if (xIsObs && meanIsObs) {
+      oocoutF(nullptr, Generation) << errPrefix << "Has two observables ?? " << std::endl;
+      return false;
+   }
+   if (!xIsObs && !meanIsObs) {
+      oocoutF(nullptr, Generation) << errPrefix << "No observable?" << std::endl;
+      return false;
    }
-   if (myobs == nullptr)  {
-      oocoutF(nullptr,Generation) << errPrefix << "No observable?" << std::endl;
+   auto *myobs = dynamic_cast<RooRealVar *>(xIsObs ? &x : &mean);
+   auto *myexp = dynamic_cast<RooAbsReal *>(xIsObs ? &mean : &x);
+   if (myobs == nullptr) {
+      oocoutF(nullptr, Generation) << errPrefix << "Observable is not a RooRealVar??" << std::endl;
       return false;
    }
    if (myexp == nullptr) {
-      oocoutF(nullptr,Generation) << errPrefix << "No observable?" << std::endl;
+      oocoutF(nullptr, Generation) << errPrefix << "Expected is not a RooAbsReal??" << std::endl;
       return false;
    }
 
@@ -989,33 +977,34 @@ bool setObsToExpected(std::span<RooAbsArg *> servers, const RooArgSet &obs, std:
 ////////////////////////////////////////////////////////////////////////////////
 /// set observed value to the expected one
 /// works for Gaussian, Poisson or LogNormal
-/// assumes mean parameter value is the argument not constant and not depending on observables
-/// (if more than two arguments are not constant will use first one but print a warning !)
 /// need to iterate on the components of the Poisson to get n and nu (nu can be a RooAbsReal)
 /// (code from G. Petrucciani and extended by L.M.)
 
-bool SetObsToExpected(RooAbsPdf &pdf, const RooArgSet &obs)
+bool SetObsToExpected(RooGaussian &pdf, const RooArgSet &obs)
 {
    std::string const &errPrefix = "AsymptoticCalculator::SetObsExpected( " + std::string{pdf.ClassName()} + " ) : ";
-   std::vector<RooAbsArg *> servers;
-   for (RooAbsArg *a : pdf.servers()) {
-      servers.emplace_back(a);
-   }
-   return setObsToExpected(servers, obs, errPrefix);
+   return setObsToExpected(const_cast<RooAbsReal &>(pdf.getX()), const_cast<RooAbsReal &>(pdf.getMean()), obs,
+                           errPrefix);
+}
+
+bool SetObsToExpected(RooPoisson &pdf, const RooArgSet &obs)
+{
+   std::string const &errPrefix = "AsymptoticCalculator::SetObsExpected( " + std::string{pdf.ClassName()} + " ) : ";
+   return setObsToExpected(const_cast<RooAbsReal &>(pdf.getX()), const_cast<RooAbsReal &>(pdf.getMean()), obs,
+                           errPrefix);
 }
 
 bool setObsToExpectedMultiVarGauss(RooMultiVarGaussian &mvgauss, const RooArgSet &obs)
 {
    // In the case of the multi-variate Gaussian, we need to iterate over the
-   // dimensions and treat the servers for each dimension separately.
+   // dimensions and treat the observable and mean for each dimension
+   // separately.
 
    std::string const &errPrefix = "AsymptoticCalculator::SetObsExpected( " + std::string{mvgauss.ClassName()} + " ) : ";
-   std::vector<RooAbsArg *> servers{nullptr, nullptr};
    bool ret = true;
    for (std::size_t iDim = 0; iDim < mvgauss.xVec().size(); ++iDim) {
-      servers[0] = &mvgauss.xVec()[iDim];
-      servers[1] = &mvgauss.muVec()[iDim];
-      ret &= setObsToExpected(servers, obs, errPrefix + " : dim " + std::to_string(iDim) + " ");
+      ret &= setObsToExpected(mvgauss.xVec()[iDim], mvgauss.muVec()[iDim], obs,
+                              errPrefix + " : dim " + std::to_string(iDim) + " ");
    }
    return ret;
 }

roofit/roostats/test/testAsymptoticCalculator.cxx

@@ -109,3 +109,41 @@ TEST(AsymptoticCalculator, SignedTestStatistic)
    std::unique_ptr<HypoTestResult> resRightSideCapped{calc.GetHypoTest()};
    EXPECT_NEAR(resRightSideCapped->NullPValue(), resRightSide->NullPValue(), 1e-6);
 }
+
+// Check that counting Asimov datasets can be generated no matter which
+// parameters are floating and even if the mean or width of a Gaussian are
+// derived quantities (covers JIRA ROOT-10096).
+TEST(AsymptoticCalculator, CountingAsimovDataSetFloatingParams)
+{
+   RooWorkspace ws;
+   ws.factory("obs[10.0, 0.0, 1000.0]");
+   ws.factory("Poisson::poisson(obs, mean[20.0, 0.0, 1000.0])");
+   ws.factory("Gaussian::gauss1(obs, mean, sigma[3.0, 1.0, 10.0])");
+   ws.factory("expr::sqrt_mean('sqrt(@0)', mean)");
+   ws.factory("Gaussian::gauss2(obs, mean, sqrt_mean)");
+   ws.factory("expr::mean2('2 * @0', mean)");
+   ws.factory("expr::sqrt_mean2('sqrt(@0)', mean2)");
+   ws.factory("Gaussian::gauss3(obs, mean2, sqrt_mean2)");
+
+   RooArgSet observables{*ws.var("obs")};
+
+   auto checkAsimov = [&](const char *pdfName, double expectedObsVal) {
+      std::unique_ptr<RooAbsData> data{
+         RooStats::AsymptoticCalculator::GenerateAsimovData(*ws.pdf(pdfName), observables)};
+      ASSERT_NE(data, nullptr) << pdfName;
+      ASSERT_EQ(data->numEntries(), 1) << pdfName;
+      EXPECT_DOUBLE_EQ(data->get(0)->getRealValue("obs"), expectedObsVal) << pdfName;
+   };
+
+   checkAsimov("poisson", 20.0);
+   // Both mean and sigma floating: used to fail with "Has two non-const arguments".
+   checkAsimov("gauss1", 20.0);
+   // Width derived from the mean: also used to fail, with no workaround for gauss3.
+   checkAsimov("gauss2", 20.0);
+   checkAsimov("gauss3", 40.0);
+
+   // With a constant mean and a floating sigma, the old server-based heuristic
+   // silently set the observable to the value of the sigma parameter.
+   ws.var("mean")->setConstant(true);
+   checkAsimov("gauss1", 20.0);
+}