Skip to content

Commit

Permalink
[RF] Fix recently-introduced RooAbsReal::getPropagatedError problems
Browse files Browse the repository at this point in the history
The recent commit `565e9aa19a` updated the `getPropagatedError` with
some additional checks, but these checks caused segfaults in some
specifit usecases:

* where the `RooAbsReal` is a parameter in the fit result itself

* where the `RooAbsReal` is a variable unrelated to the fit result

* where the `RooAbsReal` depends only on some of the parameters in the
  fit result

This commit fixes these usecases again.

A unit test that checks that these usecases don't fail anymore is now
implemented in `testRooAbsReal`.
  • Loading branch information
guitargeek committed Sep 27, 2022
1 parent 3acba96 commit df4e5f7
Show file tree
Hide file tree
Showing 2 changed files with 104 additions and 53 deletions.
48 changes: 27 additions & 21 deletions roofit/roofitcore/src/RooAbsReal.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -2605,23 +2605,36 @@ double RooAbsReal::getPropagatedError(const RooFitResult &fr, const RooArgSet &n
RooArgSet allParamsInAbsReal;
getParameters(&nset, allParamsInAbsReal);

// Strip out parameters with zero error
RooArgList paramList;
for(auto * rrvInFitResult : static_range_cast<RooRealVar*>(fr.floatParsFinal())) {
if (rrvInFitResult->getError() > 1e-20) {
auto * rrvInAbsReal = static_cast<RooRealVar*>(allParamsInAbsReal.find(*rrvInFitResult));

if(rrvInAbsReal->getVal() != rrvInFitResult->getVal()) {
throw std::runtime_error(
std::string("RooAbsReal::getPropagatedError(): the parameters of the RooAbsReal don't have") +
"the same values as in the fit result! The logic of getPropagatedError is broken in this case.");
}
for(auto * rrvFitRes : static_range_cast<RooRealVar*>(fr.floatParsFinal())) {

auto rrvInAbsReal = static_cast<RooRealVar const*>(allParamsInAbsReal.find(*rrvFitRes));

// If this RooAbsReal is a RooRealVar in the fit result, we don't need to
// propagate anything and can just return the error in the fit result
if(rrvFitRes->namePtr() == namePtr()) return rrvFitRes->getError();

// Strip out parameters with zero error
if (rrvFitRes->getError() <= 1e-20) continue;

paramList.add(*rrvInAbsReal);
// Ignore parameters in the fit result that this RooAbsReal doesn't depend on
if(!rrvInAbsReal) continue;

// Checking for float equality is a bad. We check if the values are
// negligibly far away from each other, relative to the uncertainty.
if(std::abs(rrvInAbsReal->getVal() - rrvFitRes->getVal()) > 0.01 * rrvFitRes->getError()) {
std::stringstream errMsg;
errMsg << "RooAbsReal::getPropagatedError(): the parameters of the RooAbsReal don't have"
<< " the same values as in the fit result! The logic of getPropagatedError is broken in this case.";

throw std::runtime_error(errMsg.str());
}

paramList.add(*rrvInAbsReal);
}

std::vector<double> plusVar, minusVar ;
std::vector<double> plusVar;
std::vector<double> minusVar;
plusVar.reserve(paramList.size());
minusVar.reserve(paramList.size());

Expand Down Expand Up @@ -2660,9 +2673,9 @@ double RooAbsReal::getPropagatedError(const RooFitResult &fr, const RooArgSet &n
TMatrixDSym C(paramList.getSize()) ;
vector<double> errVec(paramList.getSize()) ;
for (int i=0 ; i<paramList.getSize() ; i++) {
errVec[i] = sqrt(V(i,i)) ;
errVec[i] = std::sqrt(V(i,i)) ;
for (int j=i ; j<paramList.getSize() ; j++) {
C(i,j) = V(i,j)/sqrt(V(i,i)*V(j,j)) ;
C(i,j) = V(i,j) / std::sqrt(V(i,i)*V(j,j));
C(j,i) = C(i,j) ;
}
}
Expand All @@ -2681,13 +2694,6 @@ double RooAbsReal::getPropagatedError(const RooFitResult &fr, const RooArgSet &n










////////////////////////////////////////////////////////////////////////////////
/// Plot function or PDF on frame with support for visualization of the uncertainty encoded in the given fit result fr.
/// \param[in] frame RooPlot to plot on
Expand Down
109 changes: 77 additions & 32 deletions roofit/roofitcore/test/testRooAbsReal.cxx
Original file line number Diff line number Diff line change
@@ -1,50 +1,95 @@
// Tests for RooAbsReal
// Author: Stephan Hageboeck, CERN 05/2020
// Authors: Stephan Hageboeck, CERN 05/2020
// Jonas Rembser, CERN 09/2022

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooHelpers.h"
#include "RooGlobalFunc.h"
#include <RooAbsPdf.h>
#include <RooDataSet.h>
#include <RooFitResult.h>
#include <RooGlobalFunc.h>
#include <RooHelpers.h>
#include <RooRealVar.h>
#include <RooWorkspace.h>

#include "TTree.h"
#include "TFile.h"
#include "gtest/gtest.h"
#include <TFile.h>
#include <TTree.h>

#include <gtest/gtest.h>

#include <memory>

// ROOT-6882: Cannot read from ULong64_t branches.
TEST(RooAbsReal, ReadFromTree)
{
RooRealVar x("ULong64Branch", "xx", 0, 0, 10);
RooRealVar y("FloatBranch", "yy", 2, 0, 10);
RooRealVar z("IntBranch", "zz", 2, 0, 10);
RooRealVar a("DoubleBranch", "aa", 2, 0, 10);

TFile theFile("testRooAbsReal_1.root", "READ");
ASSERT_TRUE(theFile.IsOpen());
RooRealVar x("ULong64Branch", "xx", 0, 0, 10);
RooRealVar y("FloatBranch", "yy", 2, 0, 10);
RooRealVar z("IntBranch", "zz", 2, 0, 10);
RooRealVar a("DoubleBranch", "aa", 2, 0, 10);

TTree* tree = nullptr;
theFile.GetObject("tree", tree);
ASSERT_NE(tree, nullptr);
TFile theFile("testRooAbsReal_1.root", "READ");
ASSERT_TRUE(theFile.IsOpen());

tree->AddFriend("tree2", "testRooAbsReal_2.root");
TTree *tree = nullptr;
theFile.GetObject("tree", tree);
ASSERT_NE(tree, nullptr);

tree->AddFriend("tree2", "testRooAbsReal_2.root");

RooHelpers::HijackMessageStream hijack(RooFit::INFO, RooFit::DataHandling);
RooHelpers::HijackMessageStream hijack(RooFit::INFO, RooFit::DataHandling);

RooDataSet data("data", "data", RooArgSet(x,y,z,a), RooFit::Import(*tree));
ASSERT_EQ(data.numEntries(), 4);
RooDataSet data("data", "data", {x, y, z, a}, RooFit::Import(*tree));
ASSERT_EQ(data.numEntries(), 4);

EXPECT_DOUBLE_EQ(data.moment(x, 1), (1.+3.+3.+7.)/4.) << "ULong64Branch should contain {1, 3, 3, 7}";
EXPECT_FLOAT_EQ (data.moment(y, 1), (1.f+0.3f+0.03f+0.007f)/4.f) << "FloatBranch should contain {1., 0.3, 0.03, 0.007}";
EXPECT_DOUBLE_EQ(data.moment(z, 1), (1.+3.+3.+7.)/4.) << "IntBranch should contain {1, 3, 3, 7}";
EXPECT_DOUBLE_EQ(data.moment(a, 1), (1.+0.3+0.03+0.007)/4.) << "DoubleBranch should contain {1, 3, 3, 7}";
EXPECT_DOUBLE_EQ(data.moment(x, 1), (1. + 3. + 3. + 7.) / 4.) << "ULong64Branch should contain {1, 3, 3, 7}";
EXPECT_FLOAT_EQ(data.moment(y, 1), (1.f + 0.3f + 0.03f + 0.007f) / 4.f)
<< "FloatBranch should contain {1., 0.3, 0.03, 0.007}";
EXPECT_DOUBLE_EQ(data.moment(z, 1), (1. + 3. + 3. + 7.) / 4.) << "IntBranch should contain {1, 3, 3, 7}";
EXPECT_DOUBLE_EQ(data.moment(a, 1), (1. + 0.3 + 0.03 + 0.007) / 4.) << "DoubleBranch should contain {1, 3, 3, 7}";

const std::string& msgs = hijack.str();
const char* targetMsg = " will be converted to double precision";
EXPECT_NE(msgs.find(std::string(x.GetName()) + targetMsg), std::string::npos) << "Expect to see INFO messages for conversion of integer branch to double.";
EXPECT_NE(msgs.find(std::string(y.GetName()) + targetMsg), std::string::npos) << "Expect to see INFO messages for conversion of float branch to double.";
EXPECT_NE(msgs.find(std::string(z.GetName()) + targetMsg), std::string::npos) << "Expect to see INFO messages for conversion of integer branch to double.";
EXPECT_EQ(msgs.find(std::string(a.GetName()) + targetMsg), std::string::npos) << "Expect not to see INFO messages for conversion of double branch to double.";
const std::string &msgs = hijack.str();
const char *targetMsg = " will be converted to double precision";
EXPECT_NE(msgs.find(std::string(x.GetName()) + targetMsg), std::string::npos)
<< "Expect to see INFO messages for conversion of integer branch to double.";
EXPECT_NE(msgs.find(std::string(y.GetName()) + targetMsg), std::string::npos)
<< "Expect to see INFO messages for conversion of float branch to double.";
EXPECT_NE(msgs.find(std::string(z.GetName()) + targetMsg), std::string::npos)
<< "Expect to see INFO messages for conversion of integer branch to double.";
EXPECT_EQ(msgs.find(std::string(a.GetName()) + targetMsg), std::string::npos)
<< "Expect not to see INFO messages for conversion of double branch to double.";
}

/// Check that RooAbsReal::getPropagatedError() works as expected in some
/// corner cases.
TEST(RooAbsReal, ErrorPropagation)
{
using namespace RooFit;

RooWorkspace ws;

ws.factory("Gaussian::gauss1(x[5, 0, 10], mu[5, 0, 10], sigma1[2, 0.1, 10.])");
ws.factory("Gaussian::gauss2(x[5, 0, 10], mu[5, 0, 10], sigma2[2, 0.1, 10.])");

RooRealVar &x = *ws.var("x");
RooRealVar &mu = *ws.var("mu");
RooAbsPdf &gauss1 = *ws.pdf("gauss1");
RooAbsPdf &gauss2 = *ws.pdf("gauss2");

std::unique_ptr<RooDataSet> data{gauss1.generate(x, 1000)};

std::unique_ptr<RooFitResult> res1{gauss1.fitTo(*data, Save(), PrintLevel(-1))};

RooArgSet normSet{x};

x.setError(1.0);
EXPECT_FLOAT_EQ(x.getPropagatedError(*res1, normSet), 0.0)
<< "Propagating the uncertainties to an unrelated variable should give no uncertainty, even if that variable has "
"an intrinsic error";

EXPECT_FLOAT_EQ(mu.getPropagatedError(*res1, normSet), mu.getError())
<< "Propagating the uncertainties to an parameter in the fit result should result in that parameters uncertainty";

const double err1 = gauss1.getPropagatedError(*res1, normSet);
const double err2 = gauss2.getPropagatedError(*res1, normSet);
EXPECT_LE(err2, err1)
<< "When propagating uncertainties to another PDF that is identical to the fit model except for one parameter, "
"the uncertainty is expected to be smaller because that parameters uncertainty is not considered";
}

0 comments on commit df4e5f7

Please sign in to comment.