Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 34 additions & 45 deletions roofit/roostats/src/AsymptoticCalculator.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -919,41 +917,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;
}

Expand All @@ -969,33 +957,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;
}
Expand Down
41 changes: 41 additions & 0 deletions roofit/roostats/test/testAsymptoticCalculator.cxx
Original file line number Diff line number Diff line change
@@ -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 <memory>

// Check if asymptotic datasets for counting experiments can also be generated
// from the RooMultiVarGaussian.
TEST(AsymptoticCalculator, CountingAsimovDataSetFromMultiVarGaussian)
Expand Down Expand Up @@ -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<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);
}
Loading