diff --git a/roofit/roostats/src/AsymptoticCalculator.cxx b/roofit/roostats/src/AsymptoticCalculator.cxx index 12c9f8a611589..4af549d4b2697 100644 --- a/roofit/roostats/src/AsymptoticCalculator.cxx +++ b/roofit/roostats/src/AsymptoticCalculator.cxx @@ -69,8 +69,6 @@ The calculator can generate Asimov datasets from two kinds of PDFs: #include "TStopwatch.h" -#include - using namespace RooStats; using std::string, std::unique_ptr; @@ -919,41 +917,31 @@ void FillBins(const RooAbsPdf & pdf, const RooArgList &obs, RooAbsData & data, i } -bool setObsToExpected(std::span 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(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(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(xIsObs ? &x : &mean); + auto *myexp = dynamic_cast(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; } @@ -969,33 +957,34 @@ bool setObsToExpected(std::span 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 servers; - for (RooAbsArg *a : pdf.servers()) { - servers.emplace_back(a); - } - return setObsToExpected(servers, obs, errPrefix); + return setObsToExpected(const_cast(pdf.getX()), const_cast(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(pdf.getX()), const_cast(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 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; } diff --git a/roofit/roostats/test/testAsymptoticCalculator.cxx b/roofit/roostats/test/testAsymptoticCalculator.cxx index 6d97b7daa5c3b..abb9fe9ce7e4d 100644 --- a/roofit/roostats/test/testAsymptoticCalculator.cxx +++ b/roofit/roostats/test/testAsymptoticCalculator.cxx @@ -1,10 +1,13 @@ // Author: Jonas Rembser, CERN 01/2025 #include "RooMultiVarGaussian.h" +#include "RooRealVar.h" #include "RooStats/AsymptoticCalculator.h" #include "gtest/gtest.h" +#include + // Check if asymptotic datasets for counting experiments can also be generated // from the RooMultiVarGaussian. TEST(AsymptoticCalculator, CountingAsimovDataSetFromMultiVarGaussian) @@ -37,3 +40,41 @@ TEST(AsymptoticCalculator, CountingAsimovDataSetFromMultiVarGaussian) EXPECT_EQ(dataX.getVal(), mu.getVal()); } } + +// 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 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); +}