From 942fa47807cf43caf906b1b7ba437dd22552d05e Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sun, 6 Sep 2026 18:17:59 +0000 Subject: [PATCH 1/2] [RF] Convert stressRooStats suite to googletest MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rewrite the RooStats S.T.R.E.S.S. suite in place as a googletest suite, following the precedent of the stressRooFit conversions. All 48 test configurations of the original suite are kept, with the same models, calculator configurations, random seeds and tolerances, and the RooUnitTest semantics are reproduced (seed 12345 reset per test, failure on logged RooFit ERROR messages). Instead of comparing against the references stored in stressRooStats_ref.root, the references are now self-contained: - analytic references are computed inline in the same test (profile likelihood intervals on Gaussian/Poisson models, Bayesian central intervals via gamma quantiles) - published values from the Cousins et al. papers stay hardcoded - for the pure regression tests, the frozen values from the last stressRooStats_ref.root are hardcoded at full precision The AsymptoticCalculator test now checks both the live ProfileLikelihoodCalculator and AsymptoticCalculator significances against the frozen reference, which is strictly more coverage than the original comparison against the frozen PLC value only. The tests are parameterized over the RooFit evaluation backends via the shared gtest_wrapper.h macros, matching the backend coverage of the removed per-backend ctest invocations; on cuda builds the gtest entry takes the GPU resource lock that the dedicated cuda invocation used to declare, and both ctest entries get an explicit one-hour timeout since the backends now run serially in one process. The fixtures save and restore the global default evaluation backend so it does not leak between tests. Unlike the original suite, the default minimizer is now Minuit2 instead of Minuit. Coverage of the old minimizer is kept by a separate ctest that runs the same binary with STRESSROOSTATS_MINIMIZER=Minuit, restricted to the legacy and cpu backends like the original non-default-minimizer invocations. With the reference file and the RooUnitTest classes gone, stressRooStats_ref.root and stressRooStats_tests.h are removed (stressRooStats_models.h is still used). Also dropped: the write mode, the verbose-mode plotting code, the benchmark scaffolding, the TestHypoTestCalculator class that was never registered in any suite, and the test statistic types that were only reachable from it. 🤖 Done with the help of AI --- roofit/roostats/test/CMakeLists.txt | 23 +- roofit/roostats/test/stressRooStats.cxx | 1586 ++++++++++++---- roofit/roostats/test/stressRooStats_ref.root | Bin 29483 -> 0 bytes roofit/roostats/test/stressRooStats_tests.h | 1792 ------------------ 4 files changed, 1227 insertions(+), 2174 deletions(-) delete mode 100644 roofit/roostats/test/stressRooStats_ref.root delete mode 100644 roofit/roostats/test/stressRooStats_tests.h diff --git a/roofit/roostats/test/CMakeLists.txt b/roofit/roostats/test/CMakeLists.txt index e14c02d8a35e6..1896a906e15f9 100644 --- a/roofit/roostats/test/CMakeLists.txt +++ b/roofit/roostats/test/CMakeLists.txt @@ -6,20 +6,19 @@ ROOT_ADD_GTEST(testHypoTestInvResult testHypoTestInvResult.cxx ROOT_ADD_GTEST(testSPlot testSPlot.cxx LIBRARIES RooStats) #--stressRooStats---------------------------------------------------------------------------------- -ROOT_EXECUTABLE(stressRooStats stressRooStats.cxx LIBRARIES RooStats Gpad Net) +# Googletest version of the old RooStats S.T.R.E.S.S. suite, parameterized +# over the RooFit evaluation backends. Like in the original suite, the tests +# are also run with the Minuit2 minimizer (the default minimizer of the suite +# is Minuit2). +ROOT_ADD_GTEST(stressRooStats stressRooStats.cxx LIBRARIES RooStats LABELS longtest TIMEOUT 3600) if(mathmore) target_compile_definitions(stressRooStats PRIVATE ROOFITMORE) endif() - -configure_file(stressRooStats_ref.root stressRooStats_ref.root COPYONLY) -if(roofit_legacy_eval_backend) - ROOT_ADD_TEST(test-stressroostats-legacy COMMAND stressRooStats -b legacy FAILREGEX "FAILED|Error in" LABELS longtest) -endif() -ROOT_ADD_TEST(test-stressroostats-cpu COMMAND stressRooStats -b cpu FAILREGEX "FAILED|Error in" LABELS longtest) if(cuda) - ROOT_ADD_TEST(test-stressroostats-cuda COMMAND stressRooStats -b cuda FAILREGEX "FAILED|Error in" LABELS longtest RESOURCE_LOCK GPU) -endif() -if(roofit_legacy_eval_backend) - ROOT_ADD_TEST(test-stressroostats-legacy-minuit2 COMMAND stressRooStats -minim Minuit2 -b legacy FAILREGEX "FAILED|Error in" LABELS longtest) + set_tests_properties(gtest-roofit-roostats-stressRooStats PROPERTIES RESOURCE_LOCK GPU) endif() -ROOT_ADD_TEST(test-stressroostats-cpu-minuit2 COMMAND stressRooStats -minim Minuit2 -b cpu FAILREGEX "FAILED|Error in" LABELS longtest) +# Like the non-default-minimizer variant in the original suite, the Minuit +# variant only runs the legacy and cpu backends (that is also why it doesn't +# need to lock the GPU resource). +ROOT_ADD_TEST(test-stressroostats-minuit COMMAND stressRooStats --gtest_filter=-*EvalBackendcuda* + ENVIRONMENT STRESSROOSTATS_MINIMIZER=Minuit LABELS longtest TIMEOUT 3600) diff --git a/roofit/roostats/test/stressRooStats.cxx b/roofit/roostats/test/stressRooStats.cxx index 916fb9348eea7..f72b1bf2889a3 100644 --- a/roofit/roostats/test/stressRooStats.cxx +++ b/roofit/roostats/test/stressRooStats.cxx @@ -1,421 +1,1267 @@ -// @(#)root/roofitcore:$name: $:$id$ -// Authors: Wouter Verkerke November 2007 +// Authors: Ioan Gabriel Bucur, Lorenzo Moneta, Wouter Verkerke +// +// Googletest version of the old RooStats S.T.R.E.S.S. suite. Each of the 48 +// tests of the original suite is translated to a (parameterized) gtest test +// case with the same models, calculator configurations, random seeds and +// tolerances. +// +// The original suite compared results with references stored in +// stressRooStats_ref.root. Depending on the test, these references were either +// computed analytically, taken from a publication, or produced by the tested +// calculator itself at reference-writing time (pure regression tests). In this +// version, analytic references are computed inline, and published or +// regression reference values are hardcoded (the latter extracted from the +// last stressRooStats_ref.root). +// +// The tests are parameterized over the RooFit evaluation backends, matching +// the backend coverage of the removed stressRooStats invocations. Like the +// old "-minim" command line option, the minimizer can be chosen via the +// STRESSROOSTATS_MINIMIZER environment variable (default is Minuit2). + +#include "../../roofitcore/test/gtest_wrapper.h" + +// Global functions that build the more complex RooStats models +#include "stressRooStats_models.h" + +// RooStats headers +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include -// C/C++ headers -#include -#include -#include -#include -#include - -// Math headers -#include "Math/MinimizerOptions.h" +// RooFit headers +#include +#include +#include +#include +#include +#include +#include +#include +#include // ROOT headers -#include "TSystem.h" -#include "TString.h" -#include "TStopwatch.h" -#include "TROOT.h" -#include "TLine.h" -#include "TFile.h" -#include "TClass.h" -#include "TF1.h" -#include "TBenchmark.h" - -// RooFit headers -#include "RooGlobalFunc.h" -#include "RooNumIntConfig.h" -#include "RooMsgService.h" -#include "RooResolutionModel.h" -#include "RooRandom.h" +#include +#include +#include +#include +#include -// Tests file -#include "stressRooStats_tests.h" +#include +#include +#include +#include +#include -using std::string, std::list, std::setw, std::setfill, std::left; +using namespace ROOT::Math; using namespace RooFit; - -//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*// -// // -// RooStats Unit Test S.T.R.E.S.S. Suite // -// Authors: Ioan Gabriel Bucur, Lorenzo Moneta, Wouter Verkerke // -// // -//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*// - -//////////////////////////////////////////////////////////////////////////////// -/// Print test program number and its title - -void StatusPrint(const int id, const TString &title, const int status, const int lineWidth) +using namespace RooStats; + +namespace { + +enum ECalculatorType { kAsymptotic = 0, kFrequentist = 1, kHybrid = 2 }; +enum ETestStatType { + kSimpleLR = 0, + kRatioLR = 1, + kProfileLR = 2, + kProfileLROneSided = 3, + kProfileLROneSidedDiscovery = 4 +}; + +// Confidence levels corresponding to one, two and three Gaussian sigmas +const double kCL1Sigma = 2 * normal_cdf(1) - 1; +const double kCL2Sigma = 2 * normal_cdf(2) - 1; +const double kCL3Sigma = 2 * normal_cdf(3) - 1; + +// Value test tolerance of the original RooUnitTest +const double kVTol = 1e-3; + +std::unique_ptr buildHypoTestCalculator(const ECalculatorType calculatorType, + RooAbsData &data, const ModelConfig &nullModel, + const ModelConfig &altModel, const UInt_t toysNull, + const UInt_t toysAlt) { - TString header = TString::Format("Test %d : %s ", id, title.Data()); - std::cout << left << setw(lineWidth) << setfill('.') << header << " " - << (status > 0 ? "OK" : (status < 0 ? "SKIPPED" : "FAILED")) << std::endl; + if (calculatorType == kAsymptotic) { + return std::make_unique(data, altModel, nullModel); + } else if (calculatorType == kFrequentist) { + auto fc = std::make_unique(data, altModel, nullModel); + // set toys for speedup + fc->SetToys(toysNull, toysAlt); + return fc; + } + // kHybrid + auto hc = std::make_unique(data, altModel, nullModel); + // set toys for speedup + hc->SetToys(toysNull, toysAlt); + return hc; } -//////////////////////////////////////////////////////////////////////////////// -/// width of lines when printing test results - -int stressRooStats(const char *refFile, bool writeRef, int verbose, bool allTests, bool oneTest, int testNumber, - bool dryRun, bool doDump, bool doTreeStore) +std::unique_ptr +buildTestStatistic(const ETestStatType testStatType, const ModelConfig &nullModel, const ModelConfig &altModel) { - const int lineWidth = 120; - - // Save memory directory location - auto memDir = gDirectory; - RooUnitTest::setMemDir(gDirectory); + if (testStatType == kSimpleLR) { + auto slrts = std::make_unique(*nullModel.GetPdf(), *altModel.GetPdf()); + if (nullModel.GetSnapshot()) { + RooArgSet nullParams(*nullModel.GetSnapshot()); + if (nullModel.GetNuisanceParameters()) + nullParams.add(*nullModel.GetNuisanceParameters()); + slrts->SetNullParameters(nullParams); + } + if (altModel.GetSnapshot()) { + RooArgSet altParams(*altModel.GetSnapshot()); + if (altModel.GetNuisanceParameters()) + altParams.add(*altModel.GetNuisanceParameters()); + slrts->SetAltParameters(altParams); + } + slrts->SetAlwaysReuseNLL(true); + return slrts; + } else if (testStatType == kRatioLR) { + auto roplts = std::make_unique(*nullModel.GetPdf(), *altModel.GetPdf(), + altModel.GetSnapshot()); + roplts->SetSubtractMLE(false); + roplts->SetAlwaysReuseNLL(true); + return roplts; + } + // kProfileLR, kProfileLROneSided and kProfileLROneSidedDiscovery + auto plts = std::make_unique(*nullModel.GetPdf()); + if (testStatType == kProfileLROneSided) { + plts->SetOneSided(true); + } else if (testStatType == kProfileLROneSidedDiscovery) { + plts->SetOneSidedDiscovery(true); + } + plts->SetAlwaysReuseNLL(true); + return plts; +} - if (doTreeStore) { - RooAbsData::setDefaultStorageType(RooAbsData::Tree); +/// Create the Poisson product model in the workspace and add the observed +/// values to its data set. Returns the "S+B" model configuration (or nullptr +/// if the workspace content is unexpected, to be caught by the caller). +ModelConfig *setupPoissonProductModel(RooWorkspace &ws, int obsValueX, int obsValueY) +{ + buildPoissonProductModel(&ws); + auto model = dynamic_cast(ws.obj("S+B")); + if (model) { + ws.var("x")->setVal(obsValueX); + ws.var("y")->setVal(obsValueY); + ws.data("data")->add(*model->GetObservables()); } + return model; +} - TFile *fref = nullptr; - if (!dryRun) { - if (TString(refFile).Contains("http:")) { - if (writeRef) { - std::cout << "stressRooStats ERROR: reference file must be local file in writing mode" << std::endl; - return -1; +/// Global setup that mirrors the environment of the original stressRooStats +/// executable: minimizer selection and RooFit message streams silenced below +/// the ERROR level. +class StressRooStatsEnvironment : public ::testing::Environment { +public: + void SetUp() override + { + const char *minimizer = gSystem->Getenv("STRESSROOSTATS_MINIMIZER"); + ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizer ? minimizer : "Minuit2"); + + // Disable RooFit messages below the ERROR level, but keep a dedicated + // error stream active so that problems remain visible in the test log. + auto &msgSvc = RooMsgService::instance(); + msgSvc.setSilentMode(true); + for (int i = 0; i < msgSvc.numStreams(); ++i) { + if (msgSvc.getStream(i).minLevel < RooFit::ERROR) { + msgSvc.setStreamStatus(i, false); } - fref = TFile::Open(refFile); - } else { - fref = new TFile(refFile, writeRef ? "RECREATE" : ""); - } - if (fref->IsZombie()) { - std::cout << "stressRooStats ERROR: cannot open reference file " << refFile << std::endl; - return -1; } - } + msgSvc.addStream(RooFit::ERROR); - if (dryRun) { - // Preload singletons here so they don't show up in trace accounting - RooNumIntConfig::defaultConfig(); - RooResolutionModel::identity(); - } + AsymptoticCalculator::SetPrintLevel(0); - // Add dedicated logging stream for errors that will remain active in silent mode - RooMsgService::instance().addStream(RooFit::ERROR); - - std::cout << left << setw(lineWidth) << setfill('*') << "" << std::endl; - std::cout << "*" << setw(lineWidth - 2) << setfill(' ') << " RooStats S.T.R.E.S.S. suite " - << "*" << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - - TStopwatch timer; - timer.Start(); - - list testList; - - // 1-5 TEST PLC CONFINT SIMPLE GAUSSIAN : Confidence Level range is (0,1) - testList.push_back(new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 0.99999)); // boundary case CL -> 1 - testList.push_back( - new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 2 * ROOT::Math::normal_cdf(3) - 1)); // 3 sigma - testList.push_back( - new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 2 * ROOT::Math::normal_cdf(2) - 1)); // 2 sigma - testList.push_back( - new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 2 * ROOT::Math::normal_cdf(1) - 1)); // 1 sigma - testList.push_back(new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 0.00001)); // boundary case CL -> 0 - - // 6-10 TEST PLC CONFINT SIMPLE POISSON : Observed value range is [0,1000] - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 0)); // boundary Poisson value (0) - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 1)); - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 5)); - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 100)); - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 800)); // boundary Poisson value - - // 11-13 TEST PLC CONFINT PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for y=2*s*1.2^beta - testList.push_back(new TestProfileLikelihoodCalculator3(fref, writeRef, verbose, 10, 30)); - testList.push_back(new TestProfileLikelihoodCalculator3(fref, writeRef, verbose, 20, 25)); - testList.push_back( - new TestProfileLikelihoodCalculator3(fref, writeRef, verbose, 15, 20, 2 * ROOT::Math::normal_cdf(2) - 1)); - - // 14 TEST PLC HYPOTEST ON/OFF MODEL - testList.push_back(new TestProfileLikelihoodCalculator4(fref, writeRef, verbose)); - - // 15-18 TEST BC CONFINT CENTRAL SIMPLE POISSON : Observed value range is [0,100] - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 1)); - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 3)); - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 10)); - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 50)); - - // 19 TEST BC CONFINT SHORTEST SIMPLE POISSON - testList.push_back(new TestBayesianCalculator2(fref, writeRef, verbose)); - - // 20-22 TEST BC CONFINT CENTRAL PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for - // y=2*s*1.2^beta - testList.push_back(new TestBayesianCalculator3(fref, writeRef, verbose, 10, 30)); - testList.push_back(new TestBayesianCalculator3(fref, writeRef, verbose, 20, 25)); - testList.push_back(new TestBayesianCalculator3(fref, writeRef, verbose, 15, 20, 2 * ROOT::Math::normal_cdf(2) - 1)); - - // 23-25 TEST MCMCC CONFINT PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for y=2*s*1.2^beta - testList.push_back(new TestMCMCCalculator(fref, writeRef, verbose, 10, 30)); - testList.push_back(new TestMCMCCalculator(fref, writeRef, verbose, 20, 25)); - testList.push_back(new TestMCMCCalculator(fref, writeRef, verbose, 15, 20, 2 * ROOT::Math::normal_cdf(2) - 1)); - - // 26 TEST ZBI SIGNIFICANCE - testList.push_back(new TestZBi(fref, writeRef, verbose)); - - // 27-31 TEST PLC VS AC SIGNIFICANCE : Observed value range is [0,300] for on source and [0,1100] for off-source; tau - // has the range [0.1,5.0] - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 150, 100, 1.0)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 200, 100, 1.0)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 105, 100, 1.0)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 150, 10, 0.1)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 150, 400, 4.0)); - - // 32-36 TEST HTC SIGNIFICANCE - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kAsymptotic)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kFrequentist, kSimpleLR)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kFrequentist, kRatioLR)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kFrequentist, kProfileLROneSidedDiscovery)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kHybrid, kProfileLROneSidedDiscovery)); - - // 37-43 TEST HTI PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for y=2*s*1.2^beta - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kAsymptotic, kProfileLR, 10, 30)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kAsymptotic, kProfileLR, 20, 25)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kAsymptotic, kProfileLR, 15, 20)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kFrequentist, kProfileLR, 10, 30)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kFrequentist, kProfileLR, 20, 25)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kFrequentist, kProfileLR, 15, 20)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kHybrid, kProfileLR, 10, 30)); - - // 44-48 TEST HTI S+B+E POISSON : Observed value range is [0,50] for x = e*s+b - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kAsymptotic, kProfileLROneSided, 10, 0.95)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kAsymptotic, kProfileLROneSided, 20)); - // testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kSimpleLR, 10)); - // testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kSimpleLR, 20)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kRatioLR, 10, 0.95)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kProfileLROneSided, 10, 0.95)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kHybrid, kSimpleLR, 10, 0.95)); - - TString suiteType = TString::Format( - " Starting S.T.R.E.S.S. %s", - allTests ? "full suite" : (oneTest ? TString::Format("test %d", testNumber).Data() : "basic suite")); - - std::cout << "*" << setw(lineWidth - 3) << setfill(' ') << suiteType << " *" << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - - if (doDump) { - TFile fdbg("stressRooStats_DEBUG.root", "RECREATE"); + // NOTE: RooIntegrator1D is too slow and gives poor results +#ifdef ROOFITMORE + RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); +#endif } +}; + +[[maybe_unused]] const auto gStressRooStatsEnv = ::testing::AddGlobalTestEnvironment(new StressRooStatsEnvironment); + +// Reset random generator seeds to make results independent of test ordering, +// like RooUnitTest::runTest() did in the original suite. +void setUpStressRooStatsTest() +{ + gRandom->SetSeed(12345); + RooRandom::randomGenerator()->SetSeed(12345); + RooMsgService::instance().clearErrorCount(); +} - gBenchmark->Start("stressRooStats"); +// The original suite failed a test if RooFit ERROR messages were logged. +void tearDownStressRooStatsTest() +{ + EXPECT_EQ(RooMsgService::instance().errorCount(), 0) << "RooFit ERROR messages were logged during the test"; +} - int nFailed = 0; +/// Fixture for tests that only take the evaluation backend as parameter. The +/// RooStats calculators create their likelihoods internally, so the backend +/// parameter has to be applied via the global default backend (restored again +/// in TearDown to not leak it into other tests). +class StressRooStatsBackendTest : public ::testing::TestWithParam> { +protected: + void SetUp() override { - int i; - list::iterator iter; - - if (oneTest && (testNumber <= 0 || (UInt_t)testNumber > testList.size())) { - std::cout << "Tests are numbered from 1 to " << testList.size() << std::endl; - } else { - for (iter = testList.begin(), i = 1; iter != testList.end(); iter++, i++) { - if (!oneTest || testNumber == i) { - if (doDump) { - (*iter)->setDebug(true); - } - int status = (*iter)->isTestAvailable() ? (*iter)->runTest() : -1; - StatusPrint(i, (*iter)->GetName(), status, lineWidth); - if (!status) - nFailed++; // do not count the skipped tests - } - delete *iter; - } - } + _prevBackend = RooFit::EvalBackend::defaultValue(); + RooFit::EvalBackend::defaultValue() = std::get<0>(GetParam()).value(); + setUpStressRooStatsTest(); + } + void TearDown() override + { + tearDownStressRooStatsTest(); + RooFit::EvalBackend::defaultValue() = _prevBackend; } - gBenchmark->Stop("stressRooStats"); - - // Print table with results - bool UNIX = strcmp(gSystem->GetName(), "Unix") == 0; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - if (UNIX) { - TString sp = gSystem->GetFromPipe("uname -a"); - std::cout << "* SYS: " << sp << std::endl; - if (strstr(gSystem->GetBuildNode(), "Darwin")) { - sp = gSystem->GetFromPipe("sw_vers -productVersion"); - sp += " Mac OS X "; - std::cout << "* SYS: " << sp << std::endl; - } - } else { - const Char_t *os = gSystem->Getenv("OS"); - if (!os) { - std::cout << "* SYS: Windows 95" << std::endl; - } else { - std::cout << "* SYS: " << os << " " << gSystem->Getenv("PROCESSOR_IDENTIFIER") << std::endl; - } +private: + RooFit::EvalBackend::Value _prevBackend = RooFit::EvalBackend::Value::Cpu; +}; + +/// Fixture for tests that are parameterized over the evaluation backend and a +/// test-specific parameter set (see StressRooStatsBackendTest). +template +class StressRooStatsParamTest : public ::testing::TestWithParam> { +protected: + void SetUp() override + { + _prevBackend = RooFit::EvalBackend::defaultValue(); + RooFit::EvalBackend::defaultValue() = std::get<0>(this->GetParam()).value(); + setUpStressRooStatsTest(); + } + void TearDown() override + { + tearDownStressRooStatsTest(); + RooFit::EvalBackend::defaultValue() = _prevBackend; } - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - gBenchmark->Print("stressRooStats"); -#ifdef __CLING__ - Double_t reftime = 186.34; // pcbrun4 interpreted -#else - Double_t reftime = 93.59; // pcbrun4 compiled -#endif - const Double_t rootmarks = 860 * reftime / gBenchmark->GetCpuTime("stressRooStats"); + const ParamType ¶m() const { return std::get<1>(this->GetParam()); } + +private: + RooFit::EvalBackend::Value _prevBackend = RooFit::EvalBackend::Value::Cpu; +}; + +std::string backendName(const ::testing::TestParamInfo> &info) +{ + return "EvalBackend" + std::get<0>(info.param).name(); +} - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - std::cout << TString::Format("* ROOTMARKS = %6.1f * Root %-8s %d/%d", rootmarks, gROOT->GetVersion(), - gROOT->GetVersionDate(), gROOT->GetVersionTime()) - << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; +template +std::string backendParamName(const ::testing::TestParamInfo> &info) +{ + return "EvalBackend" + std::get<0>(info.param).name() + "_" + std::get<1>(info.param).name; +} - // NOTE: The function TStopwatch::CpuTime() calls Tstopwatch::Stop(), so you do not need to stop the timer - // separately. - std::cout << "Time at the end of job = " << timer.CpuTime() << " seconds" << std::endl; +} // namespace + +/////////////////////////////////////////////////////////////////////////////// +// +// PART ONE: PROFILE LIKELIHOOD CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - GAUSSIAN DISTRIBUTION +// +// Test the likelihood interval computed by the profile likelihood calculator +// on a Gaussian distribution. Reference interval limits are computed via +// analytic methods: solve equation 2*(ln(LL(xMax))-ln(LL(x)) = q, where q = +// normal_quantile_c(testSize/2, 1). In the case of a Gaussian distribution, +// the interval limits are equal to: +// mean +- normal_quantile_c(testSize/2, sigma/sqrt(N)). +// +// ModelConfig (implicit) : +// Observable -> x +// Parameter of Interest -> mean +// Nuisance parameter (constant!) -> sigma +// +// Original tests 1-5 (TestProfileLikelihoodCalculator1), with the confidence +// level probing the boundaries of the (0,1) range. +// +/////////////////////////////////////////////////////////////////////////////// + +struct PlcGaussianParams { + std::string name; + double confidenceLevel; +}; + +class PlcGaussianInterval : public StressRooStatsParamTest {}; + +TEST_P(PlcGaussianInterval, CompareWithAnalyticInterval) +{ + const double confidenceLevel = param().confidenceLevel; + const int N = 10; // number of observations + + // Create Gaussian model and generate a data set + RooWorkspace ws{"w"}; + ws.factory("Gaussian::gauss(x[-5,5], mean[0,-5,5], sigma[1])"); + std::unique_ptr data{ws.pdf("gauss")->generate(*ws.var("x"), N)}; + + // Reference likelihood interval limits computed via analytic methods + const double estMean = data->mean(*ws.var("x")); + const double intervalHalfWidth = + normal_quantile_c((1.0 - confidenceLevel) / 2.0, ws.var("sigma")->getValV() / std::sqrt((double)N)); + + // Calculate likelihood interval using the ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{*data, *ws.pdf("gauss"), *ws.var("mean")}; + plc.SetConfidenceLevel(confidenceLevel); + std::unique_ptr interval{plc.GetInterval()}; + + EXPECT_NEAR(interval->LowerLimit(*ws.var("mean")), estMean - intervalHalfWidth, kVTol); + EXPECT_NEAR(interval->UpperLimit(*ws.var("mean")), estMean + intervalHalfWidth, kVTol); +} - if (fref) { - fref->Close(); - delete fref; +INSTANTIATE_TEST_SUITE_P(RooStats, PlcGaussianInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PlcGaussianParams{"CLNearOne", 0.99999}, // boundary case CL -> 1 + PlcGaussianParams{"CL3Sigma", kCL3Sigma}, + PlcGaussianParams{"CL2Sigma", kCL2Sigma}, + PlcGaussianParams{"CL1Sigma", kCL1Sigma}, + PlcGaussianParams{"CLNearZero", 0.00001})), // boundary case CL -> 0 + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON DISTRIBUTION +// +// Test the 68% likelihood interval computed by the profile likelihood +// calculator on a Poisson distribution, from only one observed value. +// Reference values are computed via analytic methods: solve equation +// 2*[ln(LL(xMax)) - ln(LL(x))] = 1. +// +// ModelConfig (implicit) : +// Observable -> x +// Parameter of Interest -> mean +// +// Original tests 6-10 (TestProfileLikelihoodCalculator2), with the observed +// value probing the boundaries of the [0,1000] range. +// +/////////////////////////////////////////////////////////////////////////////// + +struct PlcPoissonParams { + std::string name; + int obsValue; +}; + +class PlcPoissonInterval : public StressRooStatsParamTest {}; + +TEST_P(PlcPoissonInterval, CompareWithAnalyticInterval) +{ + const int obsValue = param().obsValue; + + // Set a 68% confidence level for the interval + const double confidenceLevel = kCL1Sigma; + + // Create Poisson model and dataset + RooWorkspace ws{"w"}; + ws.factory(TString::Format("Poisson::poiss(x[%d,0,1000], mean[0,1000])", obsValue).Data()); + RooDataSet data{"data", "data", *ws.var("x")}; + data.add(*ws.var("x")); + + // Calculate likelihood interval using the ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{data, *ws.pdf("poiss"), *ws.var("mean")}; + plc.SetConfidenceLevel(confidenceLevel); + std::unique_ptr interval{plc.GetInterval()}; + + // Reference limits are the solutions of 2*[ln(LL(xMax)) - ln(LL(x))] = 1, + // where xMax is the point of maximum likelihood. For the special case of + // the Poisson distribution with N = 1, xMax = obsValue. + TString llRatioExpression = + TString::Format("2*(x-%d*log(x)-%d+%d*log(%d))", obsValue, obsValue, obsValue, obsValue); + // Special case obsValue = 0 because log(0) is not computable, the limit of + // n * log(n), n->0 must be taken + if (obsValue == 0) + llRatioExpression = "2*x"; + TF1 llRatio{"llRatio", llRatioExpression, 1e-100, double(obsValue)}; // lowerLimit < obsValue + + // For obsValue = 0 there is no analytic lower limit (the reference value in + // the original suite was NaN, which made its comparison pass trivially) + if (obsValue != 0) { + EXPECT_NEAR(interval->LowerLimit(*ws.var("mean")), llRatio.GetX(1), kVTol); } + llRatio.SetRange(obsValue, 1000); // upperLimit > obsValue + EXPECT_NEAR(interval->UpperLimit(*ws.var("mean")), llRatio.GetX(1), kVTol); +} - delete gBenchmark; - gBenchmark = nullptr; +INSTANTIATE_TEST_SUITE_P(RooStats, PlcPoissonInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PlcPoissonParams{"Obs0", 0}, // boundary Poisson value (0) + PlcPoissonParams{"Obs1", 1}, PlcPoissonParams{"Obs5", 5}, + PlcPoissonParams{"Obs100", 100}, + PlcPoissonParams{"Obs800", 800})), // boundary Poisson value + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON PRODUCT MODEL +// +// Test the 68% likelihood interval computed by the ProfileLikelihoodCalculator +// on a complex model. Reference values and test values are both computed with +// the ProfileLikelihoodCalculator. As such, this test can only confirm if the +// ProfileLikelihoodCalculator has the same behaviour across different computer +// platforms or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 11-13 (TestProfileLikelihoodCalculator3). +// +/////////////////////////////////////////////////////////////////////////////// + +struct PoissonProductParams { + std::string name; + int obsValueX; // observed value "x" when measuring sig + bkg1 + int obsValueY; // observed value "y" when measuring 2*sig*1.2^beta + bkg2 + double confidenceLevel; + double refLowerLimit; + double refUpperLimit; +}; + +class PlcPoissonProductInterval : public StressRooStatsParamTest {}; + +TEST_P(PlcPoissonProductInterval, RegressionInterval) +{ + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *model = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + ASSERT_NE(model, nullptr); + + // build likelihood interval with ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{*ws.data("data"), *model}; + plc.SetConfidenceLevel(param().confidenceLevel); + std::unique_ptr interval{plc.GetInterval()}; + + EXPECT_NEAR(interval->LowerLimit(*ws.var("sig")), param().refLowerLimit, kVTol); + EXPECT_NEAR(interval->UpperLimit(*ws.var("sig")), param().refUpperLimit, kVTol); +} - // Some of the object are multiple times in the list, let's make sure they - // are not deleted twice. - // The addition of memDir to the list of Cleanups is not needed if it already - // there, for example if memDir is gROOT. - bool needCleanupAdd = nullptr == gROOT->GetListOfCleanups()->FindObject(memDir->GetList()); - if (needCleanupAdd) - gROOT->GetListOfCleanups()->Add(memDir->GetList()); +INSTANTIATE_TEST_SUITE_P( + RooStats, PlcPoissonProductInterval, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PoissonProductParams{"Obs10_30_CL1Sigma", 10, 30, kCL1Sigma, 7.1428311212946145, 12.216799640536237}, + PoissonProductParams{"Obs20_25_CL1Sigma", 20, 25, kCL1Sigma, 8.7079993530346087, 14.911357163771072}, + PoissonProductParams{"Obs15_20_CL2Sigma", 15, 20, kCL2Sigma, 3.4494551349741052, 14.067570959270299})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST - ON / OFF MODEL +// +// Perform a hypothesis test using the ProfileLikelihoodCalculator on the +// on/off model. The reference values are taken from the paper: "Evaluation +// of three methods for calculating statistical significance when incorporating +// a systematic uncertainty into a test of the background-only hypothesis for +// a Poisson process" by Robert D. Cousins, James T. Linnemann, Jordan Tucker. +// +// ModelConfig (explicit) : Poisson On / Off Model +// built in stressRooStats_models.h +// +// Original test 14 (TestProfileLikelihoodCalculator4). +// +/////////////////////////////////////////////////////////////////////////////// + +class PlcOnOffHypoTest : public StressRooStatsBackendTest {}; + +TEST_P(PlcOnOffHypoTest, CompareSignificanceWithPaperValues) +{ + // A larger tolerance is needed since the values in the Cousins paper are + // given with 1e-2 precision + const double tol = 1e-2; + + // For testing purposes, we consider three special cases for which the + // values are known from the Cousins et al. paper mentioned above. The + // inputs for each of these cases are (using the notations from the paper): + // n_on, n_off, tau and Z_PL. + const int numberTestSets = 3; + const int numberOnEvents[numberTestSets] = {4, 50, 67}; + const int numberOffEvents[numberTestSets] = {5, 55, 15}; + const double tau[numberTestSets] = {5.0, 2.0, 0.5}; + const double significance[numberTestSets] = {1.95, 3.02, 3.04}; + + for (int i = 0; i < numberTestSets; ++i) { + + // build workspace and model + RooWorkspace ws{"w"}; + buildOnOffModel(ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // add observable values to data set + ws.var("n_on")->setVal(numberOnEvents[i]); + ws.var("n_off")->setVal(numberOffEvents[i]); + ws.var("tau")->setVal(tau[i]); + ws.var("tau")->setConstant(); + ws.data("data")->add(*sbModel->GetObservables()); + + // set snapshots + ws.var("sig")->setVal(numberOnEvents[i] - numberOffEvents[i] / tau[i]); + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // has as initial value a non-zero value for sig (i.e start with the S+B value) + sbModel->LoadSnapshot(); + + // get significance using the ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{*ws.data("data"), *sbModel}; + plc.SetNullParameters(*bModel->GetSnapshot()); + + std::unique_ptr result{plc.GetHypoTest()}; + EXPECT_NEAR(result->Significance(), significance[i], tol) + << "for n_on = " << numberOnEvents[i] << ", n_off = " << numberOffEvents[i] << ", tau = " << tau[i]; + } +} - memDir->GetList()->Delete("slow"); +INSTANTIATE_TEST_SUITE_P(RooStats, PlcOnOffHypoTest, ::testing::Values(ROOFIT_EVAL_BACKENDS), backendName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART TWO: BAYESIAN CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// BAYESIAN CENTRAL INTERVAL - SIMPLE POISSON MODEL +// +// Test the Bayesian central interval computed by the BayesianCalculator on a +// Poisson distribution, using different priors. The parameter of interest is +// the mean of the Poisson distribution, and there are no nuisance parameters. +// The priors used are: +// 1. constant / uniform +// 2. inverse of the mean +// 3. square root of the inverse of the mean +// 4. gamma distribution +// The posterior distribution is easily obtained analytically for these cases. +// Therefore, the reference interval limits are computed analytically. +// +// ModelConfig (implicit) : +// Observable -> x +// Parameter of Interest -> mean +// +// Original tests 15-18 (TestBayesianCalculator1). +// +/////////////////////////////////////////////////////////////////////////////// + +namespace { +double priorInvFunc(double mean) +{ + return 1.0 / mean; +} +double priorInvSqrtFunc(double mean) +{ + return 1.0 / std::sqrt(mean); +} +} // namespace + +struct BcPoissonParams { + std::string name; + int obsValue; +}; - if (needCleanupAdd) - gROOT->GetListOfCleanups()->Remove(memDir->GetList()); +class BcCentralInterval : public StressRooStatsParamTest {}; - return nFailed; +TEST_P(BcCentralInterval, CompareWithAnalyticInterval) +{ + const int obsValue = param().obsValue; + + // Set the confidence level for a 68.3% CL central interval + const double confidenceLevel = kCL1Sigma; + const double gammaShape = 2; // shape of the gamma distribution prior (gamma = alpha) + const double gammaRate = 1; // rate = 1/scale of the gamma distribution prior (beta = 1/theta) + const int numberScans = 10000; // tested to be sufficient for the scan of the Bayesian posterior + + // Create Poisson model + RooWorkspace ws{"w"}; + ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); + + // create prior pdfs + ws.factory("Uniform::prior(mean)"); + ws.import(RooCFunction1PdfBinding("priorInv", "priorInv", &priorInvFunc, *ws.var("mean"))); + ws.import( + RooCFunction1PdfBinding("priorInvSqrt", "priorInvSqrt", &priorInvSqrtFunc, *ws.var("mean"))); + ws.factory(TString::Format("Gamma::priorGamma(mean, %lf, %lf, 0)", gammaShape, gammaRate).Data()); + + // build argument sets and data set + ws.defineSet("obs", "x"); + ws.defineSet("poi", "mean"); + ws.var("x")->setVal(obsValue); + RooDataSet data{"data", "data", *ws.set("obs")}; + data.add(*ws.set("obs")); + + // Compute the interval with the BayesianCalculator for the given prior and + // compare with the analytically computed reference limits: the posterior + // for a Poisson model with the priors used here is a gamma distribution. + auto testPrior = [&](const char *priorName, double refLowerLimit, double refUpperLimit) { + BayesianCalculator bc{data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf(priorName), nullptr}; + bc.SetConfidenceLevel(confidenceLevel); + bc.SetScanOfPosterior(numberScans); + std::unique_ptr interval{bc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), refLowerLimit, kVTol) << "lower limit for prior " << priorName; + EXPECT_NEAR(interval->UpperLimit(), refUpperLimit, kVTol) << "upper limit for prior " << priorName; + }; + + const double testSize = (1.0 - confidenceLevel) / 2; + + // Uniform prior on mean + testPrior("prior", gamma_quantile(testSize, obsValue + 1, 1), // integrate to 16% + gamma_quantile_c(testSize, obsValue + 1, 1)); // integrate to 84% + // Inverse of mean prior + testPrior("priorInv", gamma_quantile(testSize, obsValue, 1), gamma_quantile_c(testSize, obsValue, 1)); + // Square root of inverse of mean prior + testPrior("priorInvSqrt", gamma_quantile(testSize, obsValue + 0.5, 1), + gamma_quantile_c(testSize, obsValue + 0.5, 1)); + // Gamma distribution prior + testPrior("priorGamma", gamma_quantile(testSize, obsValue + gammaShape, 1.0 / (1 + gammaRate)), + gamma_quantile_c(testSize, obsValue + gammaShape, 1.0 / (1 + gammaRate))); } -//_____________________________batch only_____________________ -#ifndef __CLING__ +INSTANTIATE_TEST_SUITE_P(RooStats, BcCentralInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(BcPoissonParams{"Obs1", 1}, BcPoissonParams{"Obs3", 3}, + BcPoissonParams{"Obs10", 10}, + BcPoissonParams{"Obs50", 50})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// BAYESIAN SHORTEST INTERVAL - SIMPLE POISSON MODEL +// +// Test the Bayesian shortest interval computed by the BayesianCalculator on a +// Poisson distribution, using different priors. The parameter of interest is +// the mean of the Poisson distribution, and there are no nuisance parameters. +// The priors used are: +// 1. constant / uniform +// 2. inverse of the mean +// The reference interval limits are taken from the paper: "Why isn't every +// physicist a Bayesian?" by Robert D. Cousins. +// +// Original test 19 (TestBayesianCalculator2). +// +/////////////////////////////////////////////////////////////////////////////// + +class BcShortestInterval : public StressRooStatsBackendTest {}; + +TEST_P(BcShortestInterval, CompareWithPaperValues) +{ + // the reference values in the paper have a precision of only two decimal + // points, so we increase the value tolerance accordingly + const double tol = 1e-2; + + // Put the confidence level so that we obtain a 68% confidence interval + const double confidenceLevel = kCL1Sigma; + const int obsValue = 3; // observed experiment value + const int numberScans = 100000; // sufficient number of scans + + // Create Poisson model + RooWorkspace ws{"w"}; + ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); + ws.factory("Uniform::prior(mean)"); + ws.factory("EXPR::priorInv('1/mean', mean)"); + + // build argument sets and data set + ws.defineSet("poi", "mean"); + ws.defineSet("obs", "x"); + ws.var("x")->setVal(obsValue); + RooDataSet data{"data", "data", *ws.set("obs")}; + data.add(*ws.set("obs")); + + auto testPrior = [&](const char *priorName, double refLowerLimit, double refUpperLimit) { + BayesianCalculator bc{data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf(priorName), nullptr}; + bc.SetConfidenceLevel(confidenceLevel); + bc.SetShortestInterval(); + bc.SetScanOfPosterior(numberScans); + std::unique_ptr interval{bc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), refLowerLimit, tol) << "lower limit for prior " << priorName; + EXPECT_NEAR(interval->UpperLimit(), refUpperLimit, tol) << "upper limit for prior " << priorName; + }; + + // Uniform prior on mean + testPrior("prior", 1.55, 5.15); + // Inverse of mean prior + testPrior("priorInv", 0.86, 3.85); +} -int main(int argc, const char *argv[]) +INSTANTIATE_TEST_SUITE_P(RooStats, BcShortestInterval, ::testing::Values(ROOFIT_EVAL_BACKENDS), backendName); + +/////////////////////////////////////////////////////////////////////////////// +// +// BAYESIAN CENTRAL INTERVAL - POISSON PRODUCT MODEL +// +// Test the validity of the central interval computed by the BayesianCalculator +// on a complex Poisson model distribution. Reference values and test values +// are both computed with the BayesianCalculator. As such, this test can only +// confirm if the BayesianCalculator has the same behaviour across different +// computing platforms or RooStats revisions. A uniform prior PDF is used for +// the parameter of interest ("sig"). +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 20-22 (TestBayesianCalculator3). +// +/////////////////////////////////////////////////////////////////////////////// + +class BcPoissonProductInterval : public StressRooStatsParamTest {}; + +TEST_P(BcPoissonProductInterval, RegressionInterval) { - bool doWrite = false; - int verbose = 0; - bool allTests = false; - bool oneTest = false; - int testNumber = 0; - bool dryRun = false; - bool doDump = false; - bool doTreeStore = false; - auto backend = RooFit::EvalBackend(RooFit::EvalBackend::Value::Legacy); - - // string refFileName = "http://root.cern/files/stressRooStats_v534_ref.root" ; - string refFileName = "stressRooStats_ref.root"; - string minimizerName = "Minuit"; - - // Parse command line arguments - for (int i = 1; i < argc; i++) { - string arg = argv[i]; - - if (arg == "-b") { - std::string mode = argv[++i]; - backend = RooFit::EvalBackend(mode); - std::cout << "stressRooStats: NLL evaluation backend set to " << mode << std::endl; - } else if (arg == "-f") { - std::cout << "stressRooStats: using reference file " << argv[i + 1] << std::endl; - refFileName = argv[++i]; - } else if (arg == "-w") { - std::cout << "stressRooStats: running in writing mode to update reference file" << std::endl; - doWrite = true; - } else if (arg == "-mc") { - std::cout << "stressRooStats: running in memcheck mode, no regression tests are performed" << std::endl; - dryRun = true; - } else if (arg == "-min" || arg == "-minim") { - std::cout << "stressRooStats: running using minimizer " << argv[i + 1] << std::endl; - minimizerName = argv[++i]; - } else if (arg == "-ts") { - std::cout << "stressRooStats: setting tree-based storage for datasets" << std::endl; - doTreeStore = true; - } else if (arg == "-v") { - std::cout << "stressRooStats: running in verbose mode" << std::endl; - verbose = 1; - } else if (arg == "-vv") { - std::cout << "stressRooStats: running in very verbose mode" << std::endl; - verbose = 2; - } else if (arg == "-vvv") { - std::cout << "stressRooStats: running in very very verbose mode" << std::endl; - verbose = 3; - } else if (arg == "-a") { - std::cout << "stressRooStats: deploying full suite of tests" << std::endl; - allTests = true; - } else if (arg == "-n") { - std::cout << "stressRooStats: running single test" << std::endl; - oneTest = true; - testNumber = atoi(argv[++i]); - } else if (arg == "-d") { - std::cout << "stressRooStats: setting gDebug to " << argv[i + 1] << std::endl; - gDebug = atoi(argv[++i]); - } else if (arg == "-c") { - std::cout << "stressRooStats: dumping comparison file for failed tests " << std::endl; - doDump = true; - } else if (arg == "-h" || arg == "--help") { - std::cout << R"(usage: stressRooStats [ options ] - - -b : Perform every fit in the tests with the EvalBackend() command argument, where is a string - -f : use given reference file instead of default ("stressRooStats_ref.root") - -w : write reference file, instead of reading file and running comparison tests - -n N : only run test with sequential number N - -a : run full suite of tests (default is basic suite); this overrides the -n single test option - -c : dump file stressRooStats_DEBUG.root to which results of both current result and reference for each failed test are written - -mc : memory check mode, no regression test are performed. Set this flag when running with valgrind - -min : minimizer name (default is Minuit, not Minuit2) - -vs : use vector-based storage for all datasets (default is tree-based storage) - -v/-vv : set verbose mode (show result of each regression test) or very verbose mode (show all roofit output as well) - -d N : set ROOT gDebug flag to N -)"; - return 0; - } + const int numberScans = 10; // sufficient number of scans + + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *model = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + ASSERT_NE(model, nullptr); + + // Create BayesianCalculator + BayesianCalculator bc{*ws.data("data"), *model}; + bc.SetConfidenceLevel(param().confidenceLevel); + bc.SetScanOfPosterior(numberScans); + + // Obtain confidence interval by scanning the posterior function in the + // given number of points + std::unique_ptr interval{bc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), param().refLowerLimit, kVTol); + EXPECT_NEAR(interval->UpperLimit(), param().refUpperLimit, kVTol); +} + +INSTANTIATE_TEST_SUITE_P( + RooStats, BcPoissonProductInterval, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PoissonProductParams{"Obs10_30_CL1Sigma", 10, 30, kCL1Sigma, 7.1665080051981258, 12.312785237149114}, + PoissonProductParams{"Obs20_25_CL1Sigma", 20, 25, kCL1Sigma, 8.7831170668504654, 14.874961130060584}, + PoissonProductParams{"Obs15_20_CL2Sigma", 15, 20, kCL2Sigma, 3.4603352565856857, 14.186182799724543})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART THREE: MARKOV CHAIN MONTE CARLO CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// MCMC INTERVAL CALCULATOR - POISSON PRODUCT MODEL +// +// Test the validity of the confidence interval computed by the MCMCCalculator +// on a complex Poisson model distribution. Reference values and test values +// are both computed with the MCMCCalculator. As such, this test can only +// confirm if the MCMCCalculator has the same behaviour across different +// computing platforms or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 23-25 (TestMCMCCalculator). +// +/////////////////////////////////////////////////////////////////////////////// + +class McmcPoissonProductInterval : public StressRooStatsParamTest {}; + +TEST_P(McmcPoissonProductInterval, RegressionInterval) +{ + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *model = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + ASSERT_NE(model, nullptr); + + // create and configure MCMC calculator + SequentialProposal sp{0.1}; + MCMCCalculator mcmcc{*ws.data("data"), *model}; + mcmcc.SetProposalFunction(sp); + mcmcc.SetNumIters(100000); // Metropolis-Hastings algorithm iterations + mcmcc.SetNumBurnInSteps(50); // first 50 steps to be ignored as burn-in + mcmcc.SetConfidenceLevel(param().confidenceLevel); + + // calculate the confidence interval + std::unique_ptr interval{mcmcc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(*ws.var("sig")), param().refLowerLimit, kVTol); + EXPECT_NEAR(interval->UpperLimit(*ws.var("sig")), param().refUpperLimit, kVTol); +} + +INSTANTIATE_TEST_SUITE_P( + RooStats, McmcPoissonProductInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(PoissonProductParams{"Obs10_30_CL1Sigma", 10, 30, kCL1Sigma, 7.0, 11.4}, + PoissonProductParams{"Obs20_25_CL1Sigma", 20, 25, kCL1Sigma, 9.0, + 14.600000000000001}, + PoissonProductParams{"Obs15_20_CL2Sigma", 15, 20, kCL2Sigma, + 3.4000000000000004, 13.800000000000001})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART FOUR: HYPOTHESIS TEST CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// ZBI - ON / OFF MODEL +// +// Evaluate the functionality of the top level function +// NumberCountingUtils::BinomialWithTauObsZ, which computes the significance +// of a hypothesis test via a frequentist solution. This significance, called +// ZBi, is detailed in the article "Evaluation of three methods for calculating +// statistical significance when incorporating a systematic uncertainty into a +// test of the background-only hypothesis for a Poisson process" by Robert D. +// Cousins, James T. Linnemann, Jordan Tucker. The reference values are taken +// from the paper. +// +// This computation involves no likelihood fits, so it is not parameterized +// over the evaluation backends. +// +// Original test 26 (TestZBi). +// +/////////////////////////////////////////////////////////////////////////////// + +TEST(ZBiSignificance, CompareWithPaperValues) +{ + setUpStressRooStatsTest(); + + // A larger tolerance is needed since the values in the Cousins paper are + // given with 1e-2 precision + const double tol = 1e-2; + + // For testing purposes, we consider four special cases for which the values + // are known from the Cousins et al. paper mentioned above. The inputs for + // each of these cases are (using the notations from the paper): n_on, n_off + // and tau. + const int numberTestSets = 4; + const int numberOnEvents[numberTestSets] = {4, 50, 67, 200}; + const int numberOffEvents[numberTestSets] = {5, 55, 15, 10}; + const double tau[numberTestSets] = {5.0, 2.0, 0.5, 0.1}; + const double significance[numberTestSets] = {1.66, 2.93, 2.89, 2.2}; + + for (int i = 0; i < numberTestSets; ++i) { + EXPECT_NEAR(NumberCountingUtils::BinomialWithTauObsZ(numberOnEvents[i], numberOffEvents[i], tau[i]), + significance[i], tol) + << "for n_on = " << numberOnEvents[i] << ", n_off = " << numberOffEvents[i] << ", tau = " << tau[i]; } - // if (doWrite && refFileName.find("http:") == 0) { + tearDownStressRooStatsTest(); +} - // // Locate file name part in URL and update refFileName accordingly - // char* buf = new char[refFileName.size() + 1]; - // strcpy(buf, refFileName.c_str()); - // char *ptr = strrchr(buf, '/'); - // if (!ptr) ptr = strrchr(buf, ':'); - // refFileName = ptr + 1; - // delete[] buf; +/////////////////////////////////////////////////////////////////////////////// +// +// ASYMPTOTIC CALCULATOR VS PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST +// +// This test evaluates the functionality of the AsymptoticCalculator by +// comparing the significance given from a hypothesis test on the on/off model +// with the significance given by the ProfileLikelihoodCalculator. If working +// properly, the two methods should yield identical results. On top of the +// direct comparison of the two methods, both significances are also compared +// with the frozen reference values of the original suite (which were produced +// with the ProfileLikelihoodCalculator). +// +// ModelConfig (explicit) : Poisson On / Off Model +// built in stressRooStats_models.h +// +// Original tests 27-31 (TestHypoTestCalculator1). +// +/////////////////////////////////////////////////////////////////////////////// + +struct OnOffSignificanceParams { + std::string name; + int obsValueOn; // observed value "n_on" of sig + bkg + int obsValueOff; // observed value "n_off" of tau * bkg + double tau; // parameter of the model (constant with regard to integration) + double refSignificance; +}; + +class AsymptoticVsPlcSignificance : public StressRooStatsParamTest {}; + +TEST_P(AsymptoticVsPlcSignificance, CompareSignificances) +{ + // build workspace and model, add observable values to the data set and fix + // other parameters, then make the S+B and B snapshots + auto setupWorkspace = [&](RooWorkspace &ws) -> std::pair { + buildOnOffModel(ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + if (!sbModel || !bModel) + return {nullptr, nullptr}; + + ws.var("n_on")->setVal(param().obsValueOn); + ws.var("n_off")->setVal(param().obsValueOff); + ws.var("tau")->setVal(param().tau); + ws.var("tau")->setConstant(); + ws.data("data")->add(*sbModel->GetObservables()); + ws.var("bkg")->setVal(param().obsValueOff / param().tau); + + ws.var("sig")->setVal(param().obsValueOn - param().obsValueOff / param().tau); + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0.0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + return std::make_pair(sbModel, bModel); + }; + + // Hypothesis test with the ProfileLikelihoodCalculator + double significancePlc = 0.0; + { + RooWorkspace ws{"w"}; + auto [sbModel, bModel] = setupWorkspace(ws); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + ProfileLikelihoodCalculator plc{*ws.data("data"), *sbModel}; + plc.SetNullParameters(*bModel->GetSnapshot()); + plc.SetAlternateParameters(*sbModel->GetSnapshot()); + std::unique_ptr result{plc.GetHypoTest()}; + significancePlc = result->Significance(); + } - // std::cout << "stressRooStats: WARNING running in write mode, but reference file is web file, writing local file - // instead: " - // << refFileName << std::endl; - // } + // Hypothesis test with the AsymptoticCalculator + double significanceAc = 0.0; + { + RooWorkspace ws{"w"}; + auto [sbModel, bModel] = setupWorkspace(ws); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + AsymptoticCalculator atc{*ws.data("data"), *sbModel, *bModel}; + atc.SetOneSidedDiscovery(true); + std::unique_ptr result{atc.GetHypoTest()}; + significanceAc = result->Significance(); + } - // set minimizer - ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizerName.c_str()); + EXPECT_NEAR(significancePlc, param().refSignificance, kVTol); + EXPECT_NEAR(significanceAc, param().refSignificance, kVTol); +} - // set default NLL backend - RooFit::EvalBackend::defaultValue() = backend.value(); +INSTANTIATE_TEST_SUITE_P( + RooStats, AsymptoticVsPlcSignificance, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(OnOffSignificanceParams{"Obs150_100_Tau1", 150, 100, 1.0, 3.1729725711197787}, + OnOffSignificanceParams{"Obs200_100_Tau1", 200, 100, 1.0, 5.8292198861163191}, + OnOffSignificanceParams{"Obs105_100_Tau1", 105, 100, 1.0, 0.34923144876992013}, + OnOffSignificanceParams{"Obs150_10_Tau0p1", 150, 10, 0.1, 1.3181915918803733}, + OnOffSignificanceParams{"Obs150_400_Tau4", 150, 400, 4.0, 4.0985771507906898})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// HYPOTHESIS TEST CALCULATOR TEST - SIMULTANEOUS PDF MODEL +// +// This test evaluates the functionality of the HypoTestCalculator by +// calculating the significance of the signal on a simple Simultaneous Pdf +// model with two channels. Reference values and test values are both computed +// with the HypoTestCalculator. As such, this test can only confirm if the +// HypoTestCalculator has the same behaviour across different computing +// platforms or RooStats revisions. +// +// ModelConfig (explicit) : Simultaneous Model +// built in stressRooStats_models.h +// +// Original tests 32-36 (TestHypoTestCalculator2). +// +/////////////////////////////////////////////////////////////////////////////// + +struct HypoTestCalculatorParams { + std::string name; + ECalculatorType calculatorType; + ETestStatType testStatType; + double refSignificance; +}; + +class HypoTestCalculatorSignificance : public StressRooStatsParamTest {}; + +TEST_P(HypoTestCalculatorSignificance, RegressionSignificance) +{ + // Build workspace and models + RooWorkspace ws{"w"}; + buildSimultaneousModel(&ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // set snapshots + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); // value set in model + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // the test statistic is declared before the calculator because the + // calculator's sampler will hold a raw pointer to it + std::unique_ptr testStat{buildTestStatistic(param().testStatType, *bModel, *sbModel)}; + std::unique_ptr calc{ + buildHypoTestCalculator(param().calculatorType, *ws.data("data"), *bModel, *sbModel, 500, 50)}; + if (param().calculatorType == kAsymptotic) { + static_cast(*calc).SetOneSidedDiscovery(true); + } + + // ToyMCSampler configuration + auto tmcs = static_cast(calc->GetTestStatSampler()); + tmcs->SetTestStatistic(testStat.get()); + tmcs->SetUseMultiGen(true); // speedup - gBenchmark = new TBenchmark(); - return stressRooStats(refFileName.c_str(), doWrite, verbose, allTests, oneTest, testNumber, dryRun, doDump, - doTreeStore); + std::unique_ptr result{calc->GetHypoTest()}; + EXPECT_NEAR(result->Significance(), param().refSignificance, kVTol); } -//////////////////////////////////////////////////////////////////////////////// +INSTANTIATE_TEST_SUITE_P( + RooStats, HypoTestCalculatorSignificance, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(HypoTestCalculatorParams{"Asymptotic_ProfileLROneSidedDiscovery", kAsymptotic, + kProfileLROneSidedDiscovery, 2.2499193885237001}, + HypoTestCalculatorParams{"Frequentist_SimpleLR", kFrequentist, kSimpleLR, 2.3263478740408408}, + HypoTestCalculatorParams{"Frequentist_RatioLR", kFrequentist, kRatioLR, 2.4089155458154612}, + HypoTestCalculatorParams{"Frequentist_ProfileLROneSidedDiscovery", kFrequentist, + kProfileLROneSidedDiscovery, 2.5121443279304616}, + HypoTestCalculatorParams{"Hybrid_ProfileLROneSidedDiscovery", kHybrid, + kProfileLROneSidedDiscovery, 2.2571292444862254})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART FIVE: HYPOTHESIS TEST INVERTER UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// HYPOTESTINVERTER INTERVAL - POISSON PRODUCT MODEL +// +// Test the validity of the confidence interval computed by the +// HypoTestInverter on a complex Poisson model distribution. Reference values +// and test values are both computed with the HypoTestInverter. As such, this +// test can only confirm if the HypoTestInverter has the same behaviour across +// different computing platforms or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 37-43 (TestHypoTestInverter1). +// +/////////////////////////////////////////////////////////////////////////////// + +struct HypoTestInverterParams { + std::string name; + ECalculatorType calculatorType; + ETestStatType testStatType; + int obsValueX; // observed value "x" when measuring sig + bkg1 + int obsValueY; // observed value "y" when measuring 2*sig*1.2^beta + bkg2 + double refLowerLimit; + double refUpperLimit; +}; + +class HypoTestInverterInterval : public StressRooStatsParamTest {}; + +TEST_P(HypoTestInverterInterval, RegressionInterval) +{ + const double confidenceLevel = kCL1Sigma; + + // larger value test tolerance especially when using toys (difference of + // <~ 0.1 observed between using Minuit or Minuit2) + const double tol = (param().calculatorType == kAsymptotic) ? 0.01 : 0.1; + + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *sbModel = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // set snapshots + ws.var("sig")->setVal(param().obsValueX - ws.var("bkg1")->getValV()); + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // build and configure HypoTestInverter (the test statistic is declared + // first because the calculator will hold a raw pointer to it) + std::unique_ptr testStat{buildTestStatistic(param().testStatType, *sbModel, *bModel)}; + std::unique_ptr calc{ + buildHypoTestCalculator(param().calculatorType, *ws.data("data"), *sbModel, *bModel, 100, 1)}; + HypoTestInverter hti{*calc, nullptr, 1.0 - confidenceLevel}; + hti.SetTestStatistic(*testStat); + + int nscanPoints = 10; + if (param().calculatorType == kAsymptotic) { + static_cast(*calc).SetTwoSided(); + nscanPoints = 40; + } + + hti.SetFixedScan(nscanPoints, ws.var("sig")->getMin(), ws.var("sig")->getMax()); // significant speedup + + // ToyMCSampler configuration + auto tmcs = static_cast(hti.GetHypoTestCalculator()->GetTestStatSampler()); + tmcs->SetNEventsPerToy(1); // needed because we don't have an extended pdf + tmcs->SetUseMultiGen(true); // speedup -int stressRooStats() + std::unique_ptr interval{hti.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), param().refLowerLimit, tol); + EXPECT_NEAR(interval->UpperLimit(), param().refUpperLimit, tol); +} + +INSTANTIATE_TEST_SUITE_P( + RooStats, HypoTestInverterInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(HypoTestInverterParams{"Asymptotic_ProfileLR_Obs10_30", kAsymptotic, + kProfileLR, 10, 30, 7.1386012828228216, + 12.223793911203613}, + HypoTestInverterParams{"Asymptotic_ProfileLR_Obs20_25", kAsymptotic, + kProfileLR, 20, 25, 8.7069484829792732, + 14.91407242494005}, + HypoTestInverterParams{"Asymptotic_ProfileLR_Obs15_20", kAsymptotic, + kProfileLR, 15, 20, 5.750533033549976, + 10.948013504988015}, + HypoTestInverterParams{"Frequentist_ProfileLR_Obs10_30", kFrequentist, + kProfileLR, 10, 30, 7.2115486004250906, + 12.341346144450787}, + HypoTestInverterParams{"Frequentist_ProfileLR_Obs20_25", kFrequentist, + kProfileLR, 20, 25, 8.6294497270104724, + 15.078331964885663}, + HypoTestInverterParams{"Frequentist_ProfileLR_Obs15_20", kFrequentist, + kProfileLR, 15, 20, 5.6610441994727472, + 10.87954812525739}, + HypoTestInverterParams{"Hybrid_ProfileLR_Obs10_30", kHybrid, kProfileLR, 10, + 30, 7.1079825543826782, 12.553916197221213})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// HYPOTESTINVERTER UPPER LIMIT - SIGNAL + BACKGROUND + EFFICIENCY MODEL +// +// Test the validity of the upper limit computed by the HypoTestInverter +// on a complex model distribution with signal, background and efficiency. +// Reference values and test values are both computed with the +// HypoTestInverter. As such, this test can only confirm if the +// HypoTestInverter has the same behaviour across different computing platforms +// or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Signal + Background + Efficiency +// built in stressRooStats_models.h +// +// Original tests 44-48 (TestHypoTestInverter2). +// +/////////////////////////////////////////////////////////////////////////////// + +struct HypoTestInverterUpperLimitParams { + std::string name; + ECalculatorType calculatorType; + ETestStatType testStatType; + int obsValueX; // observed value "x" when measuring sig * eff + bkg + double confidenceLevel; + double refUpperLimit; + double refExpUpperLimit; + double refExpUpperLimitMinus2; + double refExpUpperLimitMinus1; + double refExpUpperLimitPlus1; + double refExpUpperLimitPlus2; +}; + +class HypoTestInverterUpperLimit : public StressRooStatsParamTest {}; + +TEST_P(HypoTestInverterUpperLimit, RegressionUpperLimit) { - bool doWrite = false; - int verbose = 0; - bool allTests = false; - bool oneTest = false; - int testNumber = 0; - bool dryRun = false; - bool doDump = false; - bool doTreeStore = false; - string refFileName = "stressRooStats_ref.root"; - - // in interpreted mode, the minimizer is hardcoded to Minuit 1 - ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit"); - - return stressRooStats(refFileName.c_str(), doWrite, verbose, allTests, oneTest, testNumber, dryRun, doDump, - doTreeStore); + // larger value test tolerance especially when using toys (difference of + // <~ 0.1 observed between using Minuit or Minuit2) + const double tol = (param().calculatorType == kAsymptotic) ? 0.02 : 0.1; + + // Create workspace and model + RooWorkspace ws{"w"}; + buildPoissonEfficiencyModel(ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // add observed values to data set + ws.var("x")->setVal(param().obsValueX); + ws.data("data")->add(*sbModel->GetObservables()); + + // set snapshots + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // calculate upper limit with HypoTestInverter (the test statistic is + // declared first because the calculator will hold a raw pointer to it) + std::unique_ptr testStat{buildTestStatistic(param().testStatType, *sbModel, *bModel)}; + std::unique_ptr calc{ + buildHypoTestCalculator(param().calculatorType, *ws.data("data"), *sbModel, *bModel, 100, 100)}; + HypoTestInverter hti{*calc, nullptr, 1.0 - param().confidenceLevel}; + hti.SetTestStatistic(*testStat); + + int nscanPoints = 10; + if (param().calculatorType == kAsymptotic) { + static_cast(*calc).SetOneSided(true); + nscanPoints = 40; + } + + hti.SetFixedScan(nscanPoints, ws.var("sig")->getMin(), ws.var("sig")->getMax()); // significant speedup + + // needed because we have no extended pdf and the ToyMC Sampler evaluation + // returns an error + auto tmcs = static_cast(hti.GetHypoTestCalculator()->GetTestStatSampler()); + tmcs->SetNEventsPerToy(1); + tmcs->SetUseMultiGen(true); // make ToyMCSampler faster + + // calculate interval and extract observed upper limit and expected upper + // limit (+- 1, 2 sigma) + std::unique_ptr interval{hti.GetInterval()}; + EXPECT_NEAR(interval->UpperLimit(), param().refUpperLimit, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(0), param().refExpUpperLimit, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(-2), param().refExpUpperLimitMinus2, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(-1), param().refExpUpperLimitMinus1, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(1), param().refExpUpperLimitPlus1, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(2), param().refExpUpperLimitPlus2, tol); } -#endif +INSTANTIATE_TEST_SUITE_P( + RooStats, HypoTestInverterUpperLimit, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(HypoTestInverterUpperLimitParams{"Asymptotic_ProfileLROneSided_Obs10_CL0p95", + kAsymptotic, kProfileLROneSided, 10, 0.95, + 24.415175188860971, 11.71352859205691, 0.0, + 4.1446917240185863, 21.580917321964726, + 34.586630899315331}, + HypoTestInverterUpperLimitParams{"Asymptotic_ProfileLROneSided_Obs20_CL1Sigma", + kAsymptotic, kProfileLROneSided, 20, + kCL1Sigma, 35.455729668223988, + 3.8701997338910363, 0.0, 0.0, + 13.115199373487766, 24.476388830684247}, + HypoTestInverterUpperLimitParams{"Frequentist_RatioLR_Obs10_CL0p95", + kFrequentist, kRatioLR, 10, 0.95, + 27.777734839933814, 13.425920689335342, + 5.2777802875157622, 6.4814835475786641, + 20.000008863407345, 30.555555555555557}, + HypoTestInverterUpperLimitParams{"Frequentist_ProfileLROneSided_Obs10_CL0p95", + kFrequentist, kProfileLROneSided, 10, 0.95, + 23.611111111111089, 14.444453307851788, + 5.2777802875157622, 6.4814835475786641, + 20.251355120121389, 29.629641114958538}, + HypoTestInverterUpperLimitParams{"Hybrid_SimpleLR_Obs10_CL0p95", kHybrid, + kSimpleLR, 10, 0.95, 25.0, + 11.111095211258252, 5.2777802875157622, + 7.6388877462703384, 20.201996699382633, + 26.376273186118574})), + backendParamName); diff --git a/roofit/roostats/test/stressRooStats_ref.root b/roofit/roostats/test/stressRooStats_ref.root deleted file mode 100644 index ad3cf3f00c2957b9e767b66d38491bc9f676e7d4..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 29483 zcmc&-3wRVow(cZEUJ@SiekFkbA()|GGt+}WhBpE(O4Rj%40*x;Ngx>%SWpmL6@?%P z2#5%Rh=_=aELRjRpyDd7Ty<4MU3{?UMOR!=xu}=>SNCM5yQ@0WGsEQeXFwbJ{Qs${ zQ+3X%I%Q==MS$aeNaZ-LnB%JYaa`qL{QFe=B;jX$Bm7dral`J#e}Bp)&1*Yt0RLO- zZ+737l*HY&|9kwYt>+_-2%0x}!{C@bk-R#9FODCl1a zMpOI&_4nMu+5Q>D{R1Ucf&Nv2%95(82^AIjfr3Dlud<}nLH;Ti-`89!;e?;{9GB|@ zCHziWf;=ONBPv|3>eW}(YXs3NX&^onl!Qt6xFFY^a{ zWhDjWKHlLW56^iyt{0btUpwG)(&StUewyYA4yPdUqR1UMFAO?x)tM%ETc)CA>(I92 zEsZvN=1gm}EETQ29<7Ps5Ea2C%Us^MMUSSBzfcOXlxoQxm|#rR+!9}(YP?HQ+)(X@ z)D6dX^zXz{E%_cM8dJ48G?r?-Ly;ue1J(Zf>EEgr&Oggk?cYo=mTFGfS&wSG!!5f> zLmoK&?e1jR^BhaHl*vpqmTDn1wrWnnExO!X-pOT^hbG=s0kKrCrM%4qW2oj*M3>^R z9FV+2mIPUFa|ceWzy9%zoC%Oi4ai1UG0_;Rh0)lmiM+r|9xm^%`?d@}`OaFFYK^ur z!5FIXcm(SK$vXtWDR?CAz$*(zE!g+V2^LyvdnOt~HCt$k6KL}fElO$Qf!qSsvU{qQ zJ(x@bsY{G#9j)dO(qHMSNw?OUbV_+aId|822+^KW-PEB0oWdXIrj8&0tquWTcM%K$ z4@Ms(HQAAjUn!K!#^NoFLUoZ~jUw@a8-17e-gxIm`ThM&6k=vH3Z;uBis+D>ZkOB5 zd6Pf6XZ8JyKVzXZHqC!lUE;LTMRM?pjDI+9lO6K42Tqg8n`*Jf3+(A)iQ;xR6+z%# zoVR;^i_)47{{wk`-7@7R!>jR^MwHGL!%MK9>pZ+el6a?}aNZkN^>Oh}la+It3Y4Rc zV-={$I3^I?H5yQFbtHGxkK|@h21qX6jfT5*)!ld=h+QUSmhMHoH`$V6;MUYt;GSX^p4C$t(ON~fb zq?;szH&NxE;uCnE#KWbie@=Bf1X)xh7(v&SLnn+K{4r!#g$70);%`huy>4tomJ=h0 z4pDHa&5&KeO?>P1g%D0D*^I1PI_tGiQsjN2Goq4;!_B+lSaIIWX?K3rYyU~erj%^< zITO*WB#p>aQg(O*0rLmv&3*ULDHrxHgm6m9=GQYJ^`ld5SEmoCa1QHT_$7x7kH(2n z+I-EQzHjmAXvn6NZ2meE(X6CEAzCQOyWCDU2G5l});yJ6brXb#SO02m1R1WElch}8 zldAJ8`=D64ztmq;P+nBxlYLH~)1j^$>VnK#S_J5?EiD%%^?}$5oo$R#!L}N~GKp=N zz?{6#3H>6#VO|zhIB%2epAJ<|?8|~{(JnHa%k4vR@{Zub=w&aQfRlr8*PXei`(ENm zX>J?w?5nwvM7tKI^*e>U3GCqy&|c%;07nTE^z=Tm{Z&%XBo6V2rRYH0eQ zViQerh#rqqge^2HDR&H-Pvla=ByBR@8it8kH-GuvqDmf4%El`&JDzI@`f9 znMV1)1Wd1*Jutzc$O4JS@_rAz*Y%BO{}Vb;R<5p0I?b}EUL#;tvnm6&)RP65M|N>} zA3r-#>UQEXCgKibTWW}A^=u-F4yQ*F$O_SD|Gsro_K>kc*O8bH#XlOHUON_?g=ZwXT<}-m8Qu5H=5^C| z1CLtH+vOY8qo);h6^g8#%b|D#i9~l@Ry;pt>nmL$S8Gg?gbX?l?GCcg^x84dn951W zjKDcUpu;0GYitJAq&z%O=P~DNH=-=6*QMmkNOePwK1cqIhv^1!4&`ZpwtYmM3Hrr;UmyB8d&-hSQA7j`oOi_8(<@)Um4up<4R)Nu!qe-g+Mo?Qr-EUsh#~@_Yg%Uh`nfO( zk}0d`_@NPxZG$%8WD)TRj9!9zbI*SI?Vne$;5vkN4)62@lS_avM@W7NJ%$oYnWygRV%cC2|b+( z$@W3dH#c5;XVW(!eU{GhRV(cLs1ei33T*(2qC@6gFbvM?nsV;qH?8kv0d`4a0qT|1 z%+L;?2M#tMu{I4_)7E!)HL?35U0{0^nPtuyn+7qdpj^z|x*@b}@*UVUB+^L`q32GtE>I5lIlTcgPyU zAy*yb!{DM@ZL-)((nE)(q$i|gk2)Q7Dbyfhrpf+;h-M{eM5dC6WP?iPs!2^5z?71i z^KB_f4<4417_C9X<{?C?Co})WL^LZ&BQll5yc&e;7g`%;gw}={ZaA4&dlVQUy9eeK z`eR8GyMn94h3bhUa|@~A$g4(6r1mi7|1xMStlT!T#Blm>sDW%AHJ-Wo!)AkKvCw*& z_P-fu_L9UR;=CZ5x&WpFEgUr0aH_}N^C@FHX41MLu}!vzu)`hf7Isb=8#Tc{PYQKg zcp|OQCeo~qM!24VsR86t4e5P&2EfBl$iCQs2rQ98a${O+T9`L$;gVJmO#|K^8BkL& zGN99;Am$R$eI2e#Z?q!sdnVw8$bbUxw51_7EoB6G0+;u}&42sk=&3(JbU3E0;Z$$7 z*K{Lb56dQv*oDZ-AqD6`g8HWp*S2nXjHDK+7VSljHt6P1t1iU6;gchDA^zN>I7x`m zs`!Vgk92zoqF3rneHt$*(K!xe^tPj?=PEaWHg~5;IFXeAc=_8;O8Z;pE#1gNTP?OAr7^0l+h}E z3XL2FM&O_#=8iAT4Yp6Qu7Rqonc%R2M?k75WyI{o3*S2k%o-h&=_(VG|2HB6lXyD_ zhuQ=jkMcI^_>YXvH<9h_~mj!ijX)0a40I*S<%Hh$wy9C5AcAH!<@+25v@5dFZX?Yiu{J2!8t0{+7)D-tNI z#tKJRnetyjR^~-JHdf4q9kwM(Zk(g@W>>E2R+IlMtFdwhT(ZWZohp`xL$>>|M+R%l z^#9cR=zRyWED$e(X5z_MUDi>JWe=zbKg(8=CW14vKp(ODgPDE$O>-a*M9)S z>2N!Fr-U`P1NXqbk=w`t0cvFR`~8wNR-MZjX2onQoOZyO8v!G$*QvtGb~%zmv`B>W zFx(yt_$zCTw612Gp?>0KOjXg@p@>s~GG>RaA8izeOk4y4YeIGyq?L;U)B4#VNXtC5 zVx+}vE^K4Uf`n~hZ}!et&-`yIa>7pQMO|~{9%k|~2fem<(PN%%i7pXmQs7@@4jnbH z<$1EfLN%A;7ACEjnv15zY%ZKNK~_S7DLS&dWUPLloEZ%_*BWhfImmvgZjz5{blJ2; z><#l|C2YUt>(@rTL2_B)<_gm4Z}NZR8eJ4EW^*C)=EmV8OsKsMU7vYqUJiJLow~I+ z;0I&mGmXtfoI1vuXuU{BXNLk#bddx|*YcGYmRFKP37Q*2*26kRFEAOI2Uhfqn2m&U z64)#vo>kv2>rVS>Cj&g(d26*XC7AO6YUhcDlP8fKhmiBD1#!V+QC@bG-a(tf&-l=nEc;sNau;Zk~uRI$xRR(ZXA*$&Q_nY{0)zdnE`Ga zLeAD&?smdWH03|9Wpm3_1BeJdmu1tjW-q=j)&Z-!hhit0;$JIvYDj!3*~vJ$DZ3Eb zd0W-yv>xyTiL@Tm^{d);(jL6k&0!Yo`o+bLWey^?92DL|^1ja6Jr|3X5F!*k`38IR z;_6x!JxK;VY}Ub{9x}H5-Y0jF&0We`<#R@Qaat=Dr%0QmlGA}BoG4uZua|DPdFk{b zFnB~~5vsL1A2jm11lDTe#xe(F@!aH0B?mt+lg_^6#%f~6vIc>i=s>Lt=HZ+Hzxp2-w*f>bdJ6HuZP1I(SZ&ZF zOBlAgNg0wd_uPi@zx)6q>q1lAI%Kr4OemX-^u}t-gn6~sd$;vwP;4!(I75MlZ2eqy zq|>I-xnyzGY=Zb9>$us3na7m>9nW6C$ija~*`Zry(di+lo^n2JGxBe*WrN+u&}6cX zo1?NjZfv&+Y<5A5jhP^8M!^v`w4Q<)rCry`iuGWN*)Ni&<0FzC+ai_N5Rh?|{3ofxS^JCeq7@YKA7?_M}cQ>EOB>0j&f|}i4Qi2kW zy^{w{@!;jkm&=u}KxJ*H1s_4~8}>2qRwKbNku8@I##_2}%$Ac-X+s{)>s~+WNY_~h zLE%vyv4BdL^oG{SzS2x3F2aE(=mzWAK~sG&Rv~k8u=NI$ZGgzF$FSOH=rM;L?+~4cda!25 z-kaXZTDK5PD1Lbr_W0388#h0X14$c~4AU62iO>G+=R;sp9WsseT0Y*}(#VhG;}I~L zSKDB0ZLwWqYl#il1(8A`?Fona>$T*Amv}6??o^Wx?D#c)N(fsqbMQFp?(MC7K4*_G z5#X!T79#e!9>p!N&mn@5(x`gvGdfmwt&%21(!08LvLdBfW&EVvIH<4S!5P!#%IcH- zNlR(FDa00oEj83M|67m>+G2F1tSpd+w~N?RmiD9#IZ(@FMZC<&P1R4cz)Qf&<&fMi zk?ejCn4CXz$tyU%gWEN5QavSIZ7a2OmfX3)6{lW54!e+-=i~N^g0hmztIJC(+Nf7c z%*~kLFP~mgjNhsPw~T92o5ZE3aC4GYO?Z3EVQ;GqS-(HGarBUF$KGnzBWK>v8%N*L z?)}C0xY{2add(xtw(dT@a_g8OyQ)twsswNdCoNJl(yi2uluGO zhkfwE_Gf0#Ta~e^&+|J+jqH5Lo0B$m{oPM}+|#kSOq~2|-;{mr7Pr3cv!Q?dti=Ow zoO@}AanjJgjJ~%l=j*q#0w1iZ zYJYm_(W<@a^L8Fd>9O7W;ikLBOd6cUjeC4wU;n*l-CGWS&@OH47hg?$tNE|^woOkzvY_vz#kZYV^g`MjorFI&I$FAa;_~}??#aFQ^O|-3#v}S1s~tYJ z?U9uaE}QV+Q>nXp4|}oyUypn;`qS;dPyOY{{3rV6cKNtw|E_KU-%|OqTYh`P^P~Gh z>b;wnZ&hy=LFeWsshz9de}cP3)Xt6n{u9U{14TRdLQHMBR*V_-O(YN!uHvZ^vYpVR zt~1fk{l@D|SWKe6IE2L{>ia;LOiX$HdTwXY3UFCWQYPcAy*3F6p?W9kYdKg$EW!MG zfjw*|!H(-92z(#c{SM6UO<8NJizvt?0<7obCh+Geq&0~Ga=Cj!|1YXn$gVz?dvVzkMaL=rKV{WGPQmN|>EV%_BdkU}+Y+GM=d znFhODIFb<4@mo!Z-YC>k#PDtjFsZ>8TUk9*(hj3e_R8h)4K;3JD1&E;S!N_e7orA> z{0uEbjhpDp)&WaNiM|{iQmAngb(uL!NYv%y5JHWc*vi6Ll42|Oh9qj-#8Ae~5)(st zHd9OhZ*?|qQS+&uAP*kn!Ut3-XAI@k5JQc>=&PhdI+4I9=IUoi(Ux;upi%%KY>T zywzEi_Nm^pP-|H6mtw3yQ zlN3`OB#TpQ^^L4{imhys$tk9CL`aIe$`CTb9!bpAgDff0R|P^!*mI1%+7DvFW*<)# zA5)O&bTX@!jH=FK3c7%|I@4eYqN?ujN!Y82qneIID5feo2wEB17_8CuDic%v2Bd^l zwLx}*tJQ#@?dg(fww6m|6&j$Hhq@$T)n!0wz^XBzQ8hb>s|I7~mbfA>mTrkF-eT#N zxSA}MZi%b43ZW4rOMD3x;KC=W6JHv|(k_7&Pb}>cSjlA2F7X9Q2JI4Gcm%lcGFex> zlZm%F2jJ4w6G5i7(2(G%4k;mZLWU+xeEARH8oH_nkktAA#24$ZP!nH-161ljeBz61 zSg47wnE`5e;tr171eU@8uAwVj07;G9#8<8`QD4Jb-7IHSFB5x33glBGH}RDwEY!po zkN~w|>p}q7@Z}x=rK%)64N%3w!flA^4B-C1R9!G}8=#m1xD8q(0lWq;eE=p^GSC1O z4J_COs8j$fb;7>^DiBz(4NzDB*lIPyLVWCm)e8V6?8zsz3;;;eN2K s8tybc6D(1O>;X1a9fQr(F{mD|R{tlKI6Hp{j#s0;qpIH>&RXvO03j-_F#rGn diff --git a/roofit/roostats/test/stressRooStats_tests.h b/roofit/roostats/test/stressRooStats_tests.h deleted file mode 100644 index 9a958201211f9..0000000000000 --- a/roofit/roostats/test/stressRooStats_tests.h +++ /dev/null @@ -1,1792 +0,0 @@ -#include "stressRooStats_models.h" // Global functions that build complex RooStats models - -// RooStats headers -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// ROOT headers -#include -#include - -// RooFit headers -#include -#include -#include -#include -#include -#include - -// STL headers -#include - -using namespace ROOT::Math; -using namespace RooFit; -using namespace RooStats; - -// testStatType = 0 Simple Likelihood Ratio (the LEP TestStat) -// = 1 Ratio of Profiled Likelihood Ratios (the Tevatron TestStat) -// = 2 Profile Likelihood Ratio (the LHC TestStat) -// = 3 Profile Likelihood One Sided (pll = 0 if mu < mu_hat) -// = 4 Profile Likelihood Signed (pll = -pll if mu < mu_hat) -// = 5 Max Likelihood Estimate as test statistic -// = 6 Number of Observed Events as test statistic -enum ECalculatorType { kAsymptotic = 0, kFrequentist = 1, kHybrid = 2 }; -enum ETestStatType { - kSimpleLR = 0, - kRatioLR = 1, - kProfileLR = 2, - kProfileLROneSided = 3, - kProfileLROneSidedDiscovery = 4, - kProfileLRSigned = 5, - kMLE = 6, - kNObs = 7 -}; -static const char *const kECalculatorTypeString[] = {"Asymptotic", "Frequentist", "Hybrid"}; -static const char *const kETestStatTypeString[] = {"Simple-Likelihood-Ratio", - "Ratio-Of-Profiled-Likelihoods", - "Profile-Likelihood-Ratio", - "Profile-Likelihood-One-Sided", - "Profile-Likelihood-One-Sided-Discovery", - "Profile-Likelihood-Signed", - "Max-Likelihood-Estimate", - "Number-Of-Observed-Events"}; -// static const char * const kETestStatTypeString[] = { "Simple Likelihood Ratio", "Ratio Of Profiled Likelihoods", -// "Profile Likelihood Ratio", "Profile Likelihood One-Sided", "Profile Likelihood One-Sided Discovery", -// "Profile Likelihood Signed", "Max Likelihood Estimate", "Number Of Observed Events" }; -static HypoTestCalculatorGeneric *buildHypoTestCalculator(const ECalculatorType calculatorType, RooAbsData &data, - const ModelConfig &nullModel, const ModelConfig &altModel, - const UInt_t toysNull, const UInt_t toysAlt); -static TestStatistic * -buildTestStatistic(const ETestStatType testStatType, const ModelConfig &sbModel, const ModelConfig &bModel); - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART ONE: -// PROFILE LIKELIHOOD CALCULATOR UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - GAUSSIAN DISTRIBUTION -// -// Test the likelihood interval computed by the profile likelihood calculator -// on a Gaussian distribution. Reference interval limits are computed via -// analytic methods: solve equation 2*(ln(LL(xMax))-ln(LL(x)) = q, where q = -// normal_quantile_c(testSize/2, 1). In the case of a Gaussian distribution, the -// interval limits are equal to: mean +- normal_quantile_c(testSize/2, sigma/sqrt(N)). -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// Nuisance parameter (Constant !) -> sigma -// -// Input Parameters: -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 03/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator1 : public RooUnitTest { -private: - double fConfidenceLevel; - -public: - TestProfileLikelihoodCalculator1(TFile *refFile, bool writeRef, Int_t verbose, double confidenceLevel = 0.95) - : RooUnitTest("ProfileLikelihoodCalculator Interval - Gaussian Model", refFile, writeRef, verbose), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - const Int_t N = 10; // number of observations - // the compared values / objects must have the same name in write / compare modes - const TString lowerLimitString = TString::Format("tplc2_lower_limit_mean_%lf", fConfidenceLevel); - const TString upperLimitString = TString::Format("tplc2_upper_limit_mean_%lf", fConfidenceLevel); - - // TODO: see why it fails for a small number of observations - // Create Gaussian model, generate data set and define - RooWorkspace ws{"w"}; - ws.factory("Gaussian::gauss(x[-5,5], mean[0,-5,5], sigma[1])"); - std::unique_ptr data{ws.pdf("gauss")->generate(*ws.var("x"), N)}; - - if (_write == true) { - - // Calculate likelihood interval from data via analytic methods - double estMean = data->mean(*ws.var("x")); - double intervalHalfWidth = - normal_quantile_c((1.0 - fConfidenceLevel) / 2.0, ws.var("sigma")->getValV() / sqrt((double)N)); - double lowerLimit = estMean - intervalHalfWidth; - double upperLimit = estMean + intervalHalfWidth; - - // Compare the limits obtained via ProfileLikelihoodCalculator with analytically estimated values - regValue(lowerLimit, lowerLimitString); - regValue(upperLimit, upperLimitString); - - } else { - - // Calculate likelihood interval using the ProfileLikelihoodCalculator - auto plc = std::make_unique(*data, *ws.pdf("gauss"), *ws.var("mean")); - plc->SetConfidenceLevel(fConfidenceLevel); - std::unique_ptr interval{plc->GetInterval()}; - - // Register analytically computed limits in the reference file - regValue(interval->LowerLimit(*ws.var("mean")), lowerLimitString); - regValue(interval->UpperLimit(*ws.var("mean")), upperLimitString); - - plc.reset(); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON DISTRIBUTION -// -// Test the 68% likelihood interval computed by the profile likelihood calculator -// on a Poisson distribution, from only one observed value. Reference values are -// computed via analytic methods: solve equation 2*[ln(LL(xMax)) - ln(LL(x))] = 1. -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// -// Input Parameters: -// obsValue -> observed value in experiment -// -// 03/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator2 : public RooUnitTest { -private: - Int_t fObsValue; - -public: - TestProfileLikelihoodCalculator2(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValue = 5) - : RooUnitTest("ProfileLikelihoodCalculator Interval - Poisson Simple Model", refFile, writeRef, verbose), - fObsValue(obsValue){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValue < 0 || fObsValue > 1000) { - Warning("isTestAvailable", "Observed value must be in the range [0,1000]. Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // the compared values / objects must have the same name in write / compare modes - const TString lowerLimitString = TString::Format("tplc2_lower_limit_mean_%d", fObsValue); - const TString upperLimitString = TString::Format("tplc2_upper_limit_mean_%d", fObsValue); - - // write reference values - if (_write == true) { - - // Solutions of equation 2*[ln(LL(xMax)) - ln(LL(x))] = 1, where xMax is the point of maximum likelihood - // For the special case of the Poisson distribution with N = 1, xMax = obsValue - TString llRatioExpression = - TString::Format("2*(x-%d*log(x)-%d+%d*log(%d))", fObsValue, fObsValue, fObsValue, fObsValue); - // Special case fObsValue = 0 because log(0) not computable, the limit of n * log(n), n->0 must be taken - if (fObsValue == 0) - llRatioExpression = TString::Format("2*x"); - - auto llRatio = std::make_unique("llRatio", llRatioExpression, 1e-100, fObsValue); // lowerLimit < obsValue - double lowerLimit = llRatio->GetX(1); - llRatio->SetRange(fObsValue, 1000); // upperLimit > obsValue - double upperLimit = llRatio->GetX(1); - - // Compare the limits obtained via ProfileLikelihoodCalculator with the likelihood ratio analytic computations - regValue(lowerLimit, lowerLimitString); - regValue(upperLimit, upperLimitString); - - // compare with reference values - } else { - - // Set a 68% confidence level for the interval - const double confidenceLevel = 2 * normal_cdf(1) - 1.0; - - // Create Poisson model and dataset - RooWorkspace ws{"w"}; - ws.factory(TString::Format("Poisson::poiss(x[%d,0,1000], mean[0,1000])", fObsValue).Data()); - RooDataSet data{"data", "data", *ws.var("x")}; - data.add(*ws.var("x")); - - // Calculate likelihood interval using the ProfileLikelihoodCalculator - auto plc = std::make_unique(data, *ws.pdf("poiss"), *ws.var("mean")); - plc->SetConfidenceLevel(confidenceLevel); - std::unique_ptr interval{plc->GetInterval()}; - - // Register externally computed limits in the reference file - regValue(interval->LowerLimit(*ws.var("mean")), lowerLimitString); - regValue(interval->UpperLimit(*ws.var("mean")), upperLimitString); - - plc.reset(); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON PRODUCT MODEL -// -// Test the 68% likelihood interval computed by the ProfileLikelihoodCalculator -// on a complex model. Reference values and test values are both computed with -// the ProfileLikelihoodCalculator. As such, this test can only confirm if the -// ProfileLikelihoodCalculator has the same behaviour across different computer -// platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator3 : public RooUnitTest { -private: - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestProfileLikelihoodCalculator3(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueX = 15, - Int_t obsValueY = 30, double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("ProfileLikelihoodCalculator Interval - Poisson Product Model", refFile, writeRef, verbose), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Create workspace and model - RooWorkspace ws{"w"}; - buildPoissonProductModel(&ws); - auto model = dynamic_cast(ws.obj("S+B")); - - // add observed values to data set - ws.var("x")->setVal(fObsValueX); - ws.var("y")->setVal(fObsValueY); - ws.data("data")->add(*model->GetObservables()); - - std::unique_ptr initialVariables{model->GetPdf()->getVariables()}; - ws.saveSnapshot("initialVariables", *initialVariables); - - // build likelihood interval with ProfileLikelihoodCalculator - ProfileLikelihoodCalculator plc{*ws.data("data"), *model}; - plc.SetConfidenceLevel(fConfidenceLevel); - std::unique_ptr interval{plc.GetInterval()}; - - regValue(interval->LowerLimit(*ws.var("sig")), - TString::Format("tplc3_lower_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(*ws.var("sig")), - TString::Format("tplc3_upper_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - - if (_verb > 1) { - ws.loadSnapshot("initialVariables"); - ws.writeToFile(TString::Format("stressRooStats_PoissonProductModel_%d_%d.root", fObsValueX, fObsValueY)); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST - ON / OFF MODEL -// -// Perform a hypothesis test using the ProfileLikelihoodCalculator on the -// on/off model. Reference values and test values are both computed with the -// ProfileLikelihoodCalculator. As such, this test can only confirm if the -// ProfileLikelihoodCalculator has the same behaviour across different -// computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson On / Off Model -// built in stressRooStats_models.cxx -// -// For a detailed description of the on/off model, see the paper: "Evaluation -// of three methods for calculating statistical significance when incorporating -// a systematic uncertainty into a test of the background-only hypothesis for -// a Poisson process" by Robert D. Cousins, James T. Linnemann, Jordan Tucker -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator4 : public RooUnitTest { -public: - TestProfileLikelihoodCalculator4(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("ProfileLikelihoodCalculator Hypothesis Test", refFile, writeRef, verbose){}; - - // Override test value tolerance - // A larger tolerance is needed since the values in the Cousins paper are given with 1e-2 precision - double vtol() override { return 1e-2; } - - bool testCode() override - { - - // For testing purposes, we consider four special cases for which the values are known from - // the Cousins et al. paper mentioned above. The inputs for each of these cases are (using - // the notations from the paper): n_on, n_off and Z_PL. We provide a certain fixed input set - // for each case. - const Int_t numberTestSets = 3; - const Int_t numberOnEvents[numberTestSets] = {4, 50, 67}; - const Int_t numberOffEvents[numberTestSets] = {5, 55, 15}; - const double tau[numberTestSets] = {5.0, 2.0, 0.5}; - const double significance[numberTestSets] = {1.95, 3.02, 3.04}; - - for (Int_t i = 0; i < numberTestSets; ++i) { - - TString stringSignificance = - TString::Format("tplc4_significance_%d_%d_%lf", numberOnEvents[i], numberOffEvents[i], tau[i]); - - if (_write == true) { - - // register reference values from Cousins et al. paper - regValue(significance[i], stringSignificance); - - } else { - - // build workspace and model - RooWorkspace ws{"w"}; - buildOnOffModel(ws); - auto sbModel = dynamic_cast(ws.obj("S+B")); - auto bModel = dynamic_cast(ws.obj("B")); - - // add observable values to data set - ws.var("n_on")->setVal(numberOnEvents[i]); - ws.var("n_off")->setVal(numberOffEvents[i]); - ws.var("tau")->setVal(tau[i]); - ws.var("tau")->setConstant(); - ws.data("data")->add(*sbModel->GetObservables()); - - // set snapshots - ws.var("sig")->setVal(numberOnEvents[i] - numberOffEvents[i] / tau[i]); - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - ws.var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // has as initial value a non-zero value for sig (i.e start with the S+B value) - sbModel->LoadSnapshot(); - - // get significance using the ProfileLikelihoodCalculator - auto plc = std::make_unique(*ws.data("data"), *sbModel); - plc->SetNullParameters(*bModel->GetSnapshot()); - // plc->SetAlternateParameters(*sbModel->GetSnapshot()); // not needed for PLC - - regValue(plc->GetHypoTest()->Significance(), stringSignificance); - } - } - - return true; - } -}; - -// -// END OF PART ONE -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART TWO: -// BAYESIAN CALCULATOR UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// BAYESIAN CENTRAL INTERVAL - SIMPLE MODEL -// -// Test the Bayesian central interval computed by the BayesianCalculator on a -// Poisson distribution, using different priors. The parameter of interest is -// the mean of the Poisson distribution, and there are no nuisance parameters. -// The priors used are: -// 1. constant / uniform -// 2. inverse of the mean -// 3. square root of the inverse of the mean -// 4. gamma distribution -// The posterior distribution is easily obtained analytically for these cases. -// Therefore, the reference interval limits will be computed analytically. -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// -// Input Parameters: -// obsValue -> observed value in experiment -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestBayesianCalculator1 : public RooUnitTest { -private: - Int_t fObsValue; - double fConfidenceLevel; - static double priorInv(double mean) { return 1.0 / mean; } - static double priorInvSqrt(double mean) { return 1.0 / sqrt(mean); } - -public: - TestBayesianCalculator1(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValue = 3, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("BayesianCalculator Central Interval - Poisson Simple Model", refFile, writeRef, verbose), - fObsValue(obsValue), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValue < 0 || fObsValue > 100) { - Warning("isTestAvailable", "Observed value must be in the range [0,100]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Set the confidence level for a 68.3% CL central interval - const double gammaShape = 2; // shape of the gamma distribution prior (gamma = alpha) - const double gammaRate = 1; // rate = 1/scale of the gamma distribution prior (beta = 1/theta) - const Int_t numberScans = 10000; // tested to be sufficient for the scan of the Bayesian posterior - - // names of tested variables must be the same in write / comparison modes - const TString lowerLimitString = TString::Format("tbc1_lower_limit_unif_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitString = TString::Format("tbc1_upper_limit_unif_%d_%lf", fObsValue, fConfidenceLevel); - const TString lowerLimitInvString = TString::Format("tbc1_lower_limit_inv_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitInvString = TString::Format("tbc1_upper_limit_inv_%d_%lf", fObsValue, fConfidenceLevel); - const TString lowerLimitInvSqrtString = - TString::Format("tbc1_lower_limit_inv_sqrt_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitInvSqrtString = - TString::Format("tbc1_upper_limit_inv_sqrt_%d_%lf", fObsValue, fConfidenceLevel); - const TString lowerLimitGammaString = - TString::Format("tbc1_lower_limit_gamma_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitGammaString = - TString::Format("tbc1_upper_limit_gamma_%d_%lf", fObsValue, fConfidenceLevel); - - if (_write == true) { - - double lowerLimit = gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue + 1, 1); // integrate to 16% - double upperLimit = gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue + 1, 1); // integrate to 84% - double lowerLimitInv = gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue, 1); - double upperLimitInv = gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue, 1); - double lowerLimitInvSqrt = gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue + 0.5, 1); - double upperLimitInvSqrt = gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue + 0.5, 1); - double lowerLimitGamma = - gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue + gammaShape, 1.0 / (1 + gammaRate)); - double upperLimitGamma = - gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue + gammaShape, 1.0 / (1 + gammaRate)); - - // Compare the limits obtained via BayesianCalculator with quantile values - regValue(lowerLimit, lowerLimitString); - regValue(upperLimit, upperLimitString); - regValue(lowerLimitInv, lowerLimitInvString); - regValue(upperLimitInv, upperLimitInvString); - regValue(lowerLimitInvSqrt, lowerLimitInvSqrtString); - regValue(upperLimitInvSqrt, upperLimitInvSqrtString); - regValue(lowerLimitGamma, lowerLimitGammaString); - regValue(upperLimitGamma, upperLimitGammaString); - - } else { - - // Create Poisson model - RooWorkspace ws{"w"}; - ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); - // TODO: see why it does not work so well for boundary observed values {0, 100} - - // create prior pdfs - ws.factory("Uniform::prior(mean)"); - ws.import(RooCFunction1PdfBinding("priorInv", "priorInv", &priorInv, *ws.var("mean"))); - ws.import( - RooCFunction1PdfBinding("priorInvSqrt", "priorInvSqrt", priorInvSqrt, *ws.var("mean"))); - ws.factory(TString::Format("Gamma::priorGamma(mean, %lf, %lf, 0)", gammaShape, gammaRate).Data()); - - // build argument sets and data set - ws.defineSet("obs", "x"); - ws.defineSet("poi", "mean"); - ws.var("x")->setVal(fObsValue); - RooDataSet data{"data", "data", *ws.set("obs")}; - data.add(*ws.set("obs")); - - // NOTE: RooIntegrator1D is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - - std::unique_ptr bc; - std::unique_ptr interval; - - // Uniform prior on mean - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("prior"), nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitString); - regValue(interval->UpperLimit(), upperLimitString); - - // Inverse of mean prior - bc = - std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorInv"), nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitInvString); - regValue(interval->UpperLimit(), upperLimitInvString); - - // Square root of inverse of mean prior - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorInvSqrt"), - nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitInvSqrtString); - regValue(interval->UpperLimit(), upperLimitInvSqrtString); - - // Gamma distribution prior - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorGamma"), - nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitGammaString); - regValue(interval->UpperLimit(), upperLimitGammaString); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// BAYESIAN SHORTEST INTERVAL - SIMPLE POISSON MODEL -// -// Test the Bayesian shortest interval computed by the BayesianCalculator on a -// Poisson distribution, using different priors. The parameter of interest is -// the mean of the Poisson distribution, and there are no nuisance parameters. -// The priors used are: -// 1. constant / uniform -// 2. inverse of the mean -// The reference interval limits are taken from the paper: "Why isn't every -// physicist a Bayesian?" by Robert D. Cousins. -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestBayesianCalculator2 : public RooUnitTest { -public: - TestBayesianCalculator2(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("BayesianCalculator Shortest Interval - Poisson Simple Model", refFile, writeRef, verbose){}; - - // the references values in the paper have a precision of only two decimal points - // in such a situation, it is natural that we increase the value tolerance - double vtol() override { return 1e-2; } - - bool testCode() override - { - - // Put the confidence level so that we obtain a 68% confidence interval - const double confidenceLevel = 2 * normal_cdf(1) - 1; - const Int_t obsValue = 3; // observed experiment value - const Int_t numberScans = 100000; // sufficient number of scans - - // names of tested variables must be the same in write / comparison modes - const TString lowerLimitString = "tbc2_lower_limit_unif"; - const TString upperLimitString = "tbc2_upper_limit_unif"; - const TString lowerLimitInvString = "tbc2_lower_limit_inv"; - const TString upperLimitInvString = "tbc2_upper_limit_inv"; - - if (_write == true) { - - // Compare the limits obtained via BayesianCalculator with given reference values - regValue(1.55, lowerLimitString); - regValue(5.15, upperLimitString); - regValue(0.86, lowerLimitInvString); - regValue(3.85, upperLimitInvString); - - } else { - - // Create Poisson model - RooWorkspace ws{"w"}; - ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); - ws.factory("Uniform::prior(mean)"); - ws.factory("EXPR::priorInv('1/mean', mean)"); - - // build argument sets and data set - ws.defineSet("poi", "mean"); - ws.defineSet("obs", "x"); - ws.var("x")->setVal(obsValue); - RooDataSet data{"data", "data", *ws.set("obs")}; - data.add(*ws.set("obs")); - - // NOTE: RooIntegrator1D is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - // Uniform prior on mean - auto bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("prior"), nullptr); - bc->SetConfidenceLevel(confidenceLevel); - bc->SetShortestInterval(); - bc->SetScanOfPosterior(numberScans); - std::unique_ptr interval{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitString); - regValue(interval->UpperLimit(), upperLimitString); - - // Inverse of mean prior - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorInv"), nullptr); - bc->SetConfidenceLevel(confidenceLevel); - bc->SetShortestInterval(); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitInvString); - regValue(interval->UpperLimit(), upperLimitInvString); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// BAYESIAN CENTRAL INTERVAL - POISSON PRODUCT MODEL -// -// Test the validity of the central interval computed by the BayesianCalculator -// on a complex Poisson model distribution. Reference values and test values -// are both computed with the BayesianCalculator. As such, this test can only -// confirm if the BayesianCalculator has the same behaviour across different -// computing platforms or RooStats revisions. A uniform prior PDF is used for the -// parameter of interest ("sig"). -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestBayesianCalculator3 : public RooUnitTest { -private: - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestBayesianCalculator3(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueX = 15, Int_t obsValueY = 30, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("BayesianCalculator Central Interval - Poisson Product Model", refFile, writeRef, verbose), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - const Int_t numberScans = 10; // sufficient number of scans - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonProductModel(w.get()); - ModelConfig *model = static_cast(w->obj("S+B")); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->var("y")->setVal(fObsValueY); - w->data("data")->add(*model->GetObservables()); - - std::unique_ptr initialVariables{model->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // NOTE: Roo1DIntegrator is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - - // Create BayesianCalculator and - auto bc = std::make_unique(*w->data("data"), *model); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - - // Obtain confidence interval by scanning the posterior function in the given number of points - std::unique_ptr interval{bc->GetInterval()}; - regValue(interval->LowerLimit(), - TString::Format("tbc3_lower_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(), - TString::Format("tbc3_upper_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - - return true; - } -}; - -// -// END OF PART TWO -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART THREE: -// MARKOV CHAIN MONTE CARLO CALCULATOR UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// MCMC INTERVAL CALCULATOR - POISSON PRODUCT MODEL -// -// Test the validity of the confidence interval computed by the MCMCCalculator -// on a complex Poisson model distribution. Reference values and test values -// are both computed with the MCMCCalculator. As such, this test can only -// confirm if the MCMCCalculator has the same behaviour across different -// computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestMCMCCalculator : public RooUnitTest { -private: - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestMCMCCalculator(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueX = 15, Int_t obsValueY = 30, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("MCMCCalculator Interval - Poisson Product Model", refFile, writeRef, verbose), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonProductModel(w.get()); - ModelConfig *model = static_cast(w->obj("S+B")); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->var("y")->setVal(fObsValueY); - w->data("data")->add(*model->GetObservables()); - - std::unique_ptr initialVariables{model->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // NOTE: Roo1DIntegrator is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - - // create and configure MCMC calculator - auto sp = std::make_unique(0.1); - auto mcmcc = std::make_unique(*w->data("data"), *model); - mcmcc->SetProposalFunction(*sp); - mcmcc->SetNumIters(100000); // Metropolis-Hastings algorithm iterations - mcmcc->SetNumBurnInSteps(50); // first 50 steps to be ignored as burn-in - mcmcc->SetConfidenceLevel(fConfidenceLevel); - - // calculate the confidence interval - std::unique_ptr interval{mcmcc->GetInterval()}; - regValue(interval->LowerLimit(*w->var("sig")), - TString::Format("mcmcc_lower_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(*w->var("sig")), - TString::Format("mcmcc_upper_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - - return true; - } -}; - -// -// END OF PART THREE -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART FOUR: -// HYPOTHESIS TEST CALCULATOR UNIT TESTS -// - -///////////////////////////////////////////////////////////////////////// -// -// ZBI - ON / OFF MODEL -// -// Evaluate the functionality of the top level functions in RooStats -// called NumberCountingUtils::BinomialWithTauObsZ. This function -// computes the significance of a hypothesis test via a frequentist -// solution. This significance, called ZBi, is detailed in the article -// "Evaluation of three methods for calculating statistical significance -// when incorporating a systematic uncertainty into a test of the -// background-only hypothesis for a Poisson process" by Robert D. Cousins, -// James T. Linnemann, Jordan Tucker. The reference values are taken -// from the paper, as well as the On / Off model on which the test is -// evaluated. -// -// ModelConfig (implicit) : Poisson On / Off Model -// built in stressRooStats_models.cxx -// implicit in NumberCountingUtils::BinomialWithTauObsZ -// -// 05/2012 - Wouter Verkerke, Lorenzo Moneta, Ioan Gabriel Bucur -// -///////////////////////////////////////////////////////////////////////// - -class TestZBi : public RooUnitTest { -public: - TestZBi(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("ZBi Significance - On / Off Model", refFile, writeRef, verbose){}; - - // Override test value tolerance - // A larger tolerance is needed since the values in the Cousins paper are given with 1e-2 precision - double vtol() override { return 1e-2; } - - bool testCode() override - { - - // For testing purposes, we consider four special cases for which the values are known from - // the Cousins et al. paper mentioned above. The inputs for each of these cases are (using - // the notations from the paper): n_on, n_off and Z_PL. We provide a certain fixed input set - // for each case. - const Int_t numberTestSets = 4; - const Int_t numberOnEvents[numberTestSets] = {4, 50, 67, 200}; - const Int_t numberOffEvents[numberTestSets] = {5, 55, 15, 10}; - const double tau[numberTestSets] = {5.0, 2.0, 0.5, 0.1}; - const double significance[numberTestSets] = {1.66, 2.93, 2.89, 2.2}; - - for (Int_t i = 0; i < numberTestSets; ++i) { - - TString stringSignificance = - TString::Format("tzbi_significance_%d_%d_%lf", numberOnEvents[i], numberOffEvents[i], tau[i]); - - if (_write == true) { - - // register reference values from Cousins et al. paper - regValue(significance[i], stringSignificance); - - } else { - - // call top level function - regValue(NumberCountingUtils::BinomialWithTauObsZ(numberOnEvents[i], numberOffEvents[i], tau[i]), - stringSignificance); - } - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// ASYMPTOTIC CALCULATOR VS PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST -// -// This test evaluates the functionality of the AsymptoticCalculator by -// comparing the significance given from a hypothesis test on the on/off model -// with the significance given by the ProfileLikelihoodCalculator. The validity -// of the PLC hypothesis test is evaluated in TestProfileLikelihoodCalculator4. -// If working properly, the two methods should yield identical results. -// -// ModelConfig (explicit) : Poisson On / Off Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueOn -> observed value "n_on" of sig + bkg -// obsValueOff -> observed value "n_off" of tau * bkg -// tau -> parameter of the model (constant with regard to integration) -// -// 05/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestCalculator1 : public RooUnitTest { -private: - Int_t fObsValueOn; - Int_t fObsValueOff; - double fTau; - -public: - TestHypoTestCalculator1(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueOn = 150, - Int_t obsValueOff = 100, double tau = 1.0) - : RooUnitTest("AsymptoticCalculator vs ProfileLikelihoodCalculator Significance - On / Off Model", refFile, - writeRef, verbose), - fObsValueOn(obsValueOn), - fObsValueOff(obsValueOff), - fTau(tau){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueOn < 0 || fObsValueOn > 300) { - Warning("isTestAvailable", "Observed value on_source=s+b must be in the range [0,300]. Skipping test..."); - return false; - } - if (fObsValueOff < 0 || fObsValueOff > 1100) { - Warning("isTestAvailable", "Observed value off_source=tau*b must be in the range [0,1100]. Skipping test..."); - return false; - } - if (fTau < 0.1 || fTau > 5.0) { - Warning("isTestAvailable", "On/Off model parameter 'tau' must be in the range [0.1,5.0]. Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // names of tested variables must be the same in write / comparison modes - TString significanceString = TString::Format("thtc1_significance_%d_%d_%lf", fObsValueOn, fObsValueOff, fTau); - - // build workspace and model - auto w = std::make_unique("w"); - buildOnOffModel(*w); - ModelConfig *sbModel = static_cast(w->obj("S+B")); - ModelConfig *bModel = static_cast(w->obj("B")); - - // add observable values to data set and fix other parameters - w->var("n_on")->setVal(fObsValueOn); - w->var("n_off")->setVal(fObsValueOff); - w->var("tau")->setVal(fTau); - w->var("tau")->setConstant(); - w->data("data")->add(*sbModel->GetObservables()); - w->var("bkg")->setVal(fObsValueOff / fTau); - - // Make snapshots - w->var("sig")->setVal(fObsValueOn - fObsValueOff / fTau); - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - w->var("sig")->setVal(0.0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // Do hypothesis test with ProfileLikelihoodCalculator - if (_write == true) { - - auto plc = std::make_unique(*w->data("data"), *sbModel); - plc->SetNullParameters(*bModel->GetSnapshot()); - plc->SetAlternateParameters(*sbModel->GetSnapshot()); - regValue(plc->GetHypoTest()->Significance(), significanceString); - - } else { // Do hypothesis test with AsymptoticCalculator - - AsymptoticCalculator::SetPrintLevel(_verb); // disable superfluous messaging - auto atc = std::make_unique(*w->data("data"), *sbModel, *bModel); - atc->SetOneSidedDiscovery(true); - regValue(atc->GetHypoTest()->Significance(), significanceString); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// HYPOTHESIS TEST CALCULATOR TEST - SIMULTANEOUS PDF MODEL -// -// This test evaluates the functionality of the HypoTestCalculator by -// calculating the significance of the signal on a simple Simultaneous Pdf -// model with two channels. Reference values and test values are both computed -// with the HypoTestCalculator. As such, this test can only confirm if the -// HypoTestCalculator has the same behaviour across different computing -// platforms or RooStats revisions. -// -// ModelConfig (explicit) : Simultaneous Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// calculatorType -> Frequentist, Hybrid or Asymptotic -// testStatType -> Profile Likelihood Ratio, Simple Likelihood Ratio, etc... -// -// 06/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestCalculator2 : public RooUnitTest { -private: - ECalculatorType fCalculatorType; - ETestStatType fTestStatType; - -public: - TestHypoTestCalculator2(TFile *refFile, bool writeRef, Int_t verbose, ECalculatorType calculatorType = kAsymptotic, - ETestStatType testStatType = kProfileLROneSidedDiscovery) - : RooUnitTest(TString::Format("HypoTestCalculator Significance - Simultaneous Pdf - %s - %s", - kECalculatorTypeString[calculatorType], kETestStatTypeString[testStatType]), - refFile, writeRef, verbose), - fCalculatorType(calculatorType), - fTestStatType(testStatType){}; - - bool testCode() override - { - - // Build workspace and models - auto w = std::make_unique("w"); - buildSimultaneousModel(w.get()); - ModelConfig *sbModel = (ModelConfig *)w->obj("S+B"); - ModelConfig *bModel = (ModelConfig *)w->obj("B"); - - // set snapshots - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); // value set in model - w->var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - AsymptoticCalculator::SetPrintLevel(_verb); // is static (don;t care if we don't use it) - - std::unique_ptr calc{buildHypoTestCalculator(fCalculatorType, *w->data("data"), *bModel, *sbModel, 500, 50)}; - if (fCalculatorType == kAsymptotic) { - ((AsymptoticCalculator *)calc.get())->SetOneSidedDiscovery(true); - } - - // ToyMCSampler configuration - ToyMCSampler *tmcs = (ToyMCSampler *)calc->GetTestStatSampler(); - tmcs->SetTestStatistic(buildTestStatistic(fTestStatType, *bModel, *sbModel)); - tmcs->SetUseMultiGen(true); // speedup - - // Register result (test significance) - std::unique_ptr htr{calc->GetHypoTest()}; - regValue(htr->Significance(), TString::Format("thtc2_significance_%s_%s", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType])); - - // corresponding visual plots (in verbose mode) - from tutorials/roofit/roostats/StandardHypoTestDemo.C - if (_verb >= 1) { - if (fCalculatorType != kAsymptotic) { - TCanvas *c = new TCanvas("thtc2_canvas", "THTC2 Canvas"); - - c->cd(1); - HypoTestPlot *plot = new HypoTestPlot(*htr, 100); - plot->SetLogYaxis(true); - plot->Draw(); - - SamplingDistribution *altDist = htr->GetAltDistribution(); - HypoTestResult htExp("Expected result"); - htExp.Append(htr.get()); - // find quantiles in alt (S+B) distribution - double p[5]; - double q[5]; - for (Int_t i = 0; i < 5; ++i) { - double sig = -2 + i; - p[i] = ROOT::Math::normal_cdf(sig, 1); - } - std::vector values = altDist->GetSamplingDistribution(); - TMath::Quantiles(values.size(), 5, &values[0], q, p, false); - - for (Int_t i = 0; i < 5; ++i) { - htExp.SetTestStatisticData(q[i]); - double sig = -2 + i; - std::cout << "Expected p-value and significance at " << sig << " sigma = " << htExp.NullPValue() - << " significance " << htExp.Significance() << " sigma " << std::endl; - } - c->SaveAs(TString::Format("thtc2_scan_%s_%s.pdf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType])); - } else { - for (Int_t i = 0; i < 5; ++i) { - double sig = -2 + i; - double pval = - AsymptoticCalculator::GetExpectedPValues(htr->NullPValue(), htr->AlternatePValue(), -sig, false); - std::cout << "Expected p-value and significance at " << sig << " sigma = " << pval << " significance " - << ROOT::Math::normal_quantile_c(pval, 1) << " sigma " << std::endl; - } - } - } - - return true; - } -}; - -// -// END OF PART FOUR -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART FIVE: -// HYPOTHESIS TEST INVERTER UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// HYPOTESTINVERTER INTERVAL - POISSON PRODUCT MODEL -// -// Test the validity of the confidence interval computed by the HypoTestInverter -// on a complex Poisson model distribution. Reference values and test values -// are both computed with the HypoTestInverter. As such, this test can only -// confirm if the HypoTestInverter has the same behaviour across different -// computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// calculatorType -> Frequentist, Hybrid or Asymptotic -// testStatType -> Profile Likelihood Ratio, Simple Likelihood Ratio, etc... -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestInverter1 : public RooUnitTest { -private: - ECalculatorType fCalculatorType; - ETestStatType fTestStatType; - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestHypoTestInverter1(TFile *refFile, bool writeRef, Int_t verbose, ECalculatorType calculatorType = kAsymptotic, - ETestStatType testStatType = kProfileLR, Int_t obsValueX = 15, Int_t obsValueY = 30, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest(TString::Format("HypoTestInverter Interval - Poisson Product Model - %s - %s", - kECalculatorTypeString[calculatorType], kETestStatTypeString[testStatType]), - refFile, writeRef, verbose), - fCalculatorType(calculatorType), - fTestStatType(testStatType), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - // larger value test tolerance especially when using toys (difference of <~ 0.1 observed between using Minuit or - // Minuit2) - // (inherited default value is 1e-3) - double vtol() override { return (fCalculatorType == kAsymptotic) ? 0.01 : 0.1; } - - bool testCode() override - { - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonProductModel(w.get()); - ModelConfig *sbModel = (ModelConfig *)w->obj("S+B"); - ModelConfig *bModel = (ModelConfig *)w->obj("B"); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->var("y")->setVal(fObsValueY); - w->data("data")->add(*sbModel->GetObservables()); - - std::unique_ptr initialVariables{sbModel->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // set snapshots - w->var("sig")->setVal(fObsValueX - w->var("bkg1")->getValV()); - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - w->var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // build and configure HypoTestInverter - AsymptoticCalculator::SetPrintLevel(_verb); - HypoTestCalculatorGeneric *calc = - buildHypoTestCalculator(fCalculatorType, *w->data("data"), *sbModel, *bModel, 100, 1); - auto hti = std::make_unique(*calc, nullptr, 1.0 - fConfidenceLevel); - hti->SetTestStatistic(*buildTestStatistic(fTestStatType, *sbModel, *bModel)); - hti->SetVerbose(_verb); - - int nscanPoints = 10; - if (fCalculatorType == kAsymptotic) { - ((AsymptoticCalculator *)calc)->SetTwoSided(); - ((AsymptoticCalculator *)calc)->SetPrintLevel(_verb); - nscanPoints = 40; - } - - hti->SetFixedScan(nscanPoints, w->var("sig")->getMin(), w->var("sig")->getMax()); // significant speedup - - // ToyMCSampler configuration - ToyMCSampler *tmcs = (ToyMCSampler *)hti->GetHypoTestCalculator()->GetTestStatSampler(); - tmcs->SetNEventsPerToy(1); // needed because we don't have an extended pdf - tmcs->SetUseMultiGen(true); // speedup - - std::unique_ptr interval{hti->GetInterval()}; - regValue(interval->LowerLimit(), - TString::Format("thti1_lower_limit_sig_%s_%s_%d_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(), - TString::Format("thti1_upper_limit_sig_%s_%s_%d_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fObsValueY, fConfidenceLevel)); - - if (_verb >= 1) { - HypoTestInverterPlot *plot = new HypoTestInverterPlot("thti1_scan", "Two-Sided Scan", interval.get()); - TCanvas *c1 = new TCanvas("thti1_canvas", "THTI Canvas"); - c1->SetLogy(false); - plot->Draw("2CL CLB"); - c1->SaveAs(TString::Format("thti1_scan_%s_%s_%d_%d_%lf.pdf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fObsValueY, fConfidenceLevel)); - - if (_verb == 2) { - const int n = interval->ArraySize(); - if (n > 0 && interval->GetResult(0)->GetNullDistribution()) { - TCanvas *c2 = new TCanvas("thti1_teststat_dist", "HTI Test Statistic Distributions", 2); - if (n > 1) { - int ny = TMath::CeilNint(sqrt((double)n)); - int nx = TMath::CeilNint(double(n) / ny); - c2->Divide(nx, ny); - } - for (int i = 0; i < n; ++i) { - if (n > 1) - c2->cd(i + 1); - SamplingDistPlot *pl = plot->MakeTestStatPlot(i); - if (pl == nullptr) - return true; - pl->SetLogYaxis(true); - pl->Draw(); - } - c2->SaveAs(TString::Format("thti1_teststat_distrib_%s_%s_%d_%d_%lf.pdf", - kECalculatorTypeString[fCalculatorType], kETestStatTypeString[fTestStatType], - fObsValueX, fObsValueY, fConfidenceLevel)); - } - } - } - - // in case of debug write the workspace in a file - if (_verb > 1) { - w->loadSnapshot("initialVariables"); - w->writeToFile(TString::Format("stressRooStats_PoissonProductModel_%d_%d.root", fObsValueX, fObsValueY)); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// HYPOTESTINVERTER UPPER LIMIT - SIGNAL + BACKGROUND + EFFICIENCY MODEL -// -// Test the validity of the upper limit computed by the HypoTestInverter -// on a complex model distribution with signal, background and efficiency. -// Reference values and test values are both computed with the HypoTestInverter. -// As such, this test can only confirm if the HypoTestInverter has the same -// behaviour across different computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Signal + Background + Efficiency -// built in stressRooStats_models.cxx -// -/// Input Parameters: -// calculatorType -> Frequentist, Hybrid or Asymptotic -// testStatType -> Profile Likelihood Ratio, Simple Likelihood Ratio, etc... -// obsValueX -> observed value "x" when measuring sig * eff + bkg -// confidenceLevel -> Confidence Level of the upper limit we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestInverter2 : public RooUnitTest { -private: - ECalculatorType fCalculatorType; - ETestStatType fTestStatType; - Int_t fObsValueX; - double fConfidenceLevel; - -public: - TestHypoTestInverter2(TFile *refFile, bool writeRef, Int_t verbose, ECalculatorType calculatorType = kAsymptotic, - ETestStatType testStatType = kProfileLROneSided, Int_t obsValueX = 10, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest(TString::Format("HypoTestInverter Upper Limit - Poisson Efficiency Model - %s - %s", - kECalculatorTypeString[calculatorType], kETestStatTypeString[testStatType]), - refFile, writeRef, verbose), - fCalculatorType(calculatorType), - fTestStatType(testStatType), - fObsValueX(obsValueX), - fConfidenceLevel(confidenceLevel){}; - - // larger value test tolerance especially when using toys (difference of <~ 0.1 observed between using Minuit or - // Minuit2) - // (inherited default value is 1e-3) - double vtol() override { return (fCalculatorType == kAsymptotic) ? 0.02 : 0.1; } - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 50) { - Warning("isTestAvailable", "Observed value X=s*e+b must be in the range [0,70]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonEfficiencyModel(*w); - ModelConfig *sbModel = (ModelConfig *)w->obj("S+B"); - ModelConfig *bModel = (ModelConfig *)w->obj("B"); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->data("data")->add(*sbModel->GetObservables()); - - std::unique_ptr initialVariables{sbModel->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // set snapshots - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - w->var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // calculate upper limit with HypoTestInverter - AsymptoticCalculator::SetPrintLevel(_verb); - HypoTestCalculatorGeneric *calc = - buildHypoTestCalculator(fCalculatorType, *w->data("data"), *sbModel, *bModel, 100, 100); - auto hti = std::make_unique(*calc, nullptr, 1.0 - fConfidenceLevel); - hti->SetTestStatistic(*buildTestStatistic(fTestStatType, *sbModel, *bModel)); - hti->SetVerbose(_verb); - - int npoints = 10; - if (fCalculatorType == kAsymptotic) { - ((AsymptoticCalculator *)calc)->SetOneSided(true); - ((AsymptoticCalculator *)calc)->SetPrintLevel(_verb); - npoints = 40; - } - - hti->SetFixedScan(npoints, w->var("sig")->getMin(), w->var("sig")->getMax()); // significant speedup - - // needed because we have no extended pdf and the ToyMC Sampler evaluation returns an error - ToyMCSampler *tmcs = (ToyMCSampler *)hti->GetHypoTestCalculator()->GetTestStatSampler(); - tmcs->SetNEventsPerToy(1); - tmcs->SetUseMultiGen(true); // make ToyMCSampler faster - - // calculate interval and extract observed upper limit and expected upper limit (+- sigma) - std::unique_ptr interval{hti->GetInterval()}; - regValue(interval->UpperLimit(), - TString::Format("thti2_upper_limit_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(0), - TString::Format("thti2_exp_upper_limit_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(-2), - TString::Format("thti2_exp_upper_limit_-2_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(-1), - TString::Format("thti2_exp_upper_limit_-1_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(1), - TString::Format("thti2_exp_upper_limit_+1_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(2), - TString::Format("thti2_exp_upper_limit_+2_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - - if (_verb >= 1) { - HypoTestInverterPlot *plot = new HypoTestInverterPlot("thti2_scan", "HTI Upper Limit Scan", interval.get()); - TCanvas *c1 = new TCanvas("HypoTestInverter Scan"); - c1->SetLogy(false); - plot->Draw("2CL CLB"); - c1->SaveAs(TString::Format("thti2_scan_%s_%s_%d_%lf.pdf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - - if (_verb == 2) { - const int n = interval->ArraySize(); - if (n > 0 && interval->GetResult(0)->GetNullDistribution()) { - TCanvas *c2 = new TCanvas("thti2_teststat_dist", "HTI Test Statistic Distributions", 2); - if (n > 1) { - int ny = TMath::CeilNint(sqrt((double)n)); - int nx = TMath::CeilNint(double(n) / ny); - c2->Divide(nx, ny); - } - for (int i = 0; i < n; ++i) { - if (n > 1) - c2->cd(i + 1); - SamplingDistPlot *pl = plot->MakeTestStatPlot(i); - if (pl == nullptr) - return true; - pl->SetLogYaxis(true); - pl->Draw(); - } - c2->SaveAs(TString::Format("thti2_teststat_distrib_%s_%s_%d_%lf.pdf", - kECalculatorTypeString[fCalculatorType], kETestStatTypeString[fTestStatType], - fObsValueX, fConfidenceLevel)); - } - } - } - - if (_verb > 1) { - w->loadSnapshot("initialVariables"); - w->writeToFile("stressRooStats_PoissonEfficiencyModel.root"); - } - - return true; - } -}; - -// -// END OF PART FIVE -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -// Other tests currently not included in any suite - -class TestHypoTestCalculator : public RooUnitTest { -public: - TestHypoTestCalculator(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("HypoTestCalculator - On / Off Problem", refFile, writeRef, verbose){}; - - bool testCode() override - { - - const Int_t xValue = 150; - const Int_t yValue = 100; - const double tauValue = 1.0; - - if (_write == true) { - - // register analytical Z_Bi value - double Z_Bi = NumberCountingUtils::BinomialWithTauObsZ(xValue, yValue, tauValue); - regValue(Z_Bi, "thtc_significance_hybrid"); - - } else { - - // Make model for prototype on/off problem - // Pois(x | s+b) * Pois(y | tau b ) - auto w = std::make_unique("w"); - w->factory( - TString::Format("Poisson::on_pdf(x[%d,0,500],sum::splusb(sig[0,0,100],bkg[100,0,300]))", xValue).Data()); - w->factory(TString::Format("Poisson::off_pdf(y[%d,0,500],prod::taub(tau[%lf],bkg))", yValue, tauValue).Data()); - w->factory("PROD::prod_pdf(on_pdf, off_pdf)"); - - w->var("x")->setVal(xValue); - w->var("y")->setVal(yValue); - w->var("y")->setConstant(); - w->var("tau")->setVal(tauValue); - - // construct the Bayesian-averaged model (eg. a projection pdf) - // p'(x|s) = \int db p(x|s+b) * [ p(y|b) * prior(b) ] - w->factory("Uniform::prior(bkg)"); - w->factory("PROJ::averagedModel(PROD::foo(on_pdf|bkg,off_pdf,prior),bkg)"); - - // define sets of variables obs={x} and poi={sig} - // x is the only observable in the main measurement and y is treated as a separate measurement, - // which is used to produce the prior that will be used in the calculation to randomize the nuisance parameters - w->defineSet("obs", "x"); - w->defineSet("poi", "sig"); - - // Add observable value to a data set - RooDataSet data{"data", "data", *w->set("obs")}; - data.add(*w->set("obs")); - - // Build S+B and B models - ModelConfig *sbModel = new ModelConfig("SB_ModelConfig", w.get()); - sbModel->SetPdf(*w->pdf("prod_pdf")); - sbModel->SetObservables(*w->set("obs")); - sbModel->SetParametersOfInterest(*w->set("poi")); - w->var("sig")->setVal(xValue - yValue / tauValue); // important ! - sbModel->SetSnapshot(*w->set("poi")); - - ModelConfig *bModel = new ModelConfig("B_ModelConfig", w.get()); - bModel->SetPdf(*w->pdf("prod_pdf")); - bModel->SetObservables(*w->set("obs")); - bModel->SetParametersOfInterest(*w->set("poi")); - w->var("sig")->setVal(0.0); // important ! - bModel->SetSnapshot(*w->set("poi")); - - // alternate priors - w->factory("Gaussian::gauss_prior(bkg, y, expr::sqrty('sqrt(y)', y))"); - w->factory("Lognormal::lognorm_prior(bkg, y, expr::kappa('1+1./sqrt(y)',y))"); - - // build test statistic - SimpleLikelihoodRatioTestStat *slrts = - new SimpleLikelihoodRatioTestStat(*bModel->GetPdf(), *sbModel->GetPdf()); - slrts->SetNullParameters(*bModel->GetSnapshot()); - slrts->SetAltParameters(*sbModel->GetSnapshot()); - slrts->SetAlwaysReuseNLL(true); - - RatioOfProfiledLikelihoodsTestStat *roplts = - new RatioOfProfiledLikelihoodsTestStat(*bModel->GetPdf(), *sbModel->GetPdf()); - roplts->SetAlwaysReuseNLL(true); - - ProfileLikelihoodTestStat *pllts = new ProfileLikelihoodTestStat(*bModel->GetPdf()); - pllts->SetAlwaysReuseNLL(true); - - MaxLikelihoodEstimateTestStat *mlets = new MaxLikelihoodEstimateTestStat( - *sbModel->GetPdf(), *((RooRealVar *)sbModel->GetParametersOfInterest()->first())); - - NumEventsTestStat *nevts = new NumEventsTestStat(*sbModel->GetPdf()); - - auto htc = std::make_unique(data, *sbModel, *bModel); - ToyMCSampler *tmcs = (ToyMCSampler *)htc->GetTestStatSampler(); - tmcs->SetNEventsPerToy(1); - htc->SetToys(5000, 1000); - htc->ForcePriorNuisanceAlt(*w->pdf("off_pdf")); - htc->ForcePriorNuisanceNull(*w->pdf("off_pdf")); - - tmcs->SetTestStatistic(pllts); - std::unique_ptr htr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "PLLTS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(mlets); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "MLETS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(nevts); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "NEVTS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(slrts); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "SLRTS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(roplts); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "ROPLTS " << htr->Significance() << std::endl; - - regValue(htr->Significance(), "thtc_significance_hybrid"); - - if (_verb > 1) - w->writeToFile("stressRooStats_OnOffModel.root"); - } - - return true; - } -}; - -static HypoTestCalculatorGeneric *buildHypoTestCalculator(const ECalculatorType calculatorType, RooAbsData &data, - const ModelConfig &nullModel, const ModelConfig &altModel, - const UInt_t toysNull, const UInt_t toysAlt) -{ - HypoTestCalculatorGeneric *calc = nullptr; - - if (calculatorType == kAsymptotic) { - AsymptoticCalculator *ac = new AsymptoticCalculator(data, altModel, nullModel); - calc = ac; - } else if (calculatorType == kFrequentist) { - FrequentistCalculator *fc = new FrequentistCalculator(data, altModel, nullModel); - // set toys for speedup - fc->SetToys(toysNull, toysAlt); - calc = fc; - } else { // kHybrid - HybridCalculator *hc = new HybridCalculator(data, altModel, nullModel); - // set toys for speedup - hc->SetToys(toysNull, toysAlt); - calc = hc; - } - - assert(calc != NULL); // sanity check - should never happen - - return calc; -} - -static TestStatistic * -buildTestStatistic(const ETestStatType testStatType, const ModelConfig &nullModel, const ModelConfig &altModel) -{ - - TestStatistic *testStat = nullptr; - - if (testStatType == kSimpleLR) { - auto *slrts = new SimpleLikelihoodRatioTestStat(*nullModel.GetPdf(), *altModel.GetPdf()); - // TODO - different for HypoTestInverter and HypoTestCalculator - RooArgSet nullParams(*nullModel.GetSnapshot()); - if (nullModel.GetNuisanceParameters()) - nullParams.add(*nullModel.GetNuisanceParameters()); - if (nullModel.GetSnapshot()) - slrts->SetNullParameters(nullParams); - RooArgSet altParams(*altModel.GetSnapshot()); - if (altModel.GetNuisanceParameters()) - altParams.add(*altModel.GetNuisanceParameters()); - if (altModel.GetSnapshot()) - slrts->SetAltParameters(altParams); - slrts->SetAlwaysReuseNLL(true); - testStat = slrts; - } else if (testStatType == kRatioLR) { - auto *roplts = - new RatioOfProfiledLikelihoodsTestStat(*nullModel.GetPdf(), *altModel.GetPdf(), altModel.GetSnapshot()); - roplts->SetSubtractMLE(false); - roplts->SetAlwaysReuseNLL(true); - testStat = roplts; - } else if (testStatType == kMLE) { - auto *mlets = new MaxLikelihoodEstimateTestStat(*nullModel.GetPdf(), - *((RooRealVar *)nullModel.GetParametersOfInterest()->first())); - testStat = mlets; - } else if (testStatType == kNObs) { - NumEventsTestStat *nevtts = new NumEventsTestStat(*nullModel.GetPdf()); - testStat = nevtts; - } else { // kProfileLR, kProfileLROneSided and kProfileLRSigned - auto *plts = new ProfileLikelihoodTestStat(*nullModel.GetPdf()); - if (testStatType == kProfileLROneSided) { - plts->SetOneSided(true); - } else if (testStatType == kProfileLROneSidedDiscovery) { - plts->SetOneSidedDiscovery(true); - } else if (testStatType == kProfileLRSigned) { - plts->SetSigned(true); - } - plts->SetAlwaysReuseNLL(true); - testStat = plts; - } - - assert(testStat != nullptr); // sanity check - should never happen - - return testStat; -} From 293d3b1cd380f63359c72d3a57bd93e463ee62c3 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sun, 6 Sep 2026 19:05:16 +0000 Subject: [PATCH 2/2] [RF] Avoid RooUnitTest dependency in TestStatistics/testPlot.cxx MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Convert the RooFit::TestStatistics plotting test to a reference-free gtest, following the precedent of the other stress test conversions. The model, parallel minimization setup and likelihood plotting are unchanged, but instead of comparing the RooPlot against the reference in TestStatistics_ref.root via RooUnitTest::regPlot, the test now checks: - the fitted parameter against the analytic maximum likelihood estimate (with a tolerance at the scale of the Minuit convergence criterion) - every point of the plotted curve against a direct evaluation of the likelihood at that parameter value - the likelihood shape against the analytically known Poisson -log L (up to the constant offset) The TestStatistics_ref.root file had no other users, so it is removed together with the COPY_TO_BUILDDIR of the ctest registration. 🤖 Done with the help of AI --- roofit/roofitcore/test/CMakeLists.txt | 3 +- .../TestStatistics/TestStatistics_ref.root | Bin 12076 -> 0 bytes .../test/TestStatistics/testPlot.cxx | 142 ++++++++++-------- 3 files changed, 80 insertions(+), 65 deletions(-) delete mode 100644 roofit/roofitcore/test/TestStatistics/TestStatistics_ref.root diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 810eb310a7fb4..34d6338740196 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -94,8 +94,7 @@ ROOT_ADD_GTEST(testRooTruthModel testRooTruthModel.cxx LIBRARIES RooFitCore RooF COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/rooAbsAnaConvPdf_classV3.root) if (roofit_multiprocess) - ROOT_ADD_GTEST(testTestStatisticsPlot TestStatistics/testPlot.cxx LIBRARIES RooFitMultiProcess RooFitCore RooFit - COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/TestStatistics/TestStatistics_ref.root) + ROOT_ADD_GTEST(testTestStatisticsPlot TestStatistics/testPlot.cxx LIBRARIES RooFitMultiProcess RooFitCore RooFit) ROOT_ADD_GTEST(testLikelihoodGradientJob TestStatistics/testLikelihoodGradientJob.cxx LIBRARIES RooFitMultiProcess RooFitCore m ROOT::TestSupport) target_include_directories(testLikelihoodGradientJob PRIVATE ${RooFitCore_MultiProcess_TestStatistics_INCLUDE_DIR}) ROOT_ADD_GTEST(testLikelihoodJob TestStatistics/testLikelihoodJob.cxx LIBRARIES RooFitMultiProcess RooFitCore m) diff --git a/roofit/roofitcore/test/TestStatistics/TestStatistics_ref.root b/roofit/roofitcore/test/TestStatistics/TestStatistics_ref.root deleted file mode 100644 index 6b8c89029b9eacf7a3b97c54f31bc83c52970c01..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12076 zcmb`Mbx>SE*X9RzC%C)2Gq}6E1a}`axVyW%yF+kyf&@))cPE4(AIZ11wSVotTlIC- z>3gdC&ePAS>AL;9GcJyft^mM2B>(_00|1yAKCaInq4E*nAE68R(U}7P2z3AexgP*5 zjXP@CivaDsd5AQ`BC+@CWA;C4X@I{%I(8P5SO5V4*hl`c8UO&LE@^FNj<0SGa#eFR zcC`k%TAPACyO>*ib};u~a{24+|2heP`m1PXG0pcOw*3*hAMsBa0DzzUf7(^~=L+-R z=l|Q)-`dUo16O}b0F(QO<-Zv-g9m2{;QjsG^&|TJi9&velKgk6U`qd~=*N8@P3oTj z{Qn`4@G<%yyIT3{<)dZ)lac&aoW?)DD}QC+)si@Nd zGXBq+>|Z5tmk%K`7*$6{WjjaLPk;BrN3LpaY$yMB;Iq?5hbBVuzvf9CEl9}h-Nm%pM$tb?(QPMLckLN@5X3kW;8MrOboIMR!L0y@AOB6dLRz;P15(h&D+`9uulF06&B|UGs(7{b!5=b@wFm{+Hq!tO+P0nif@JSOD zj87;1f%*Mbp(HlUB1DPb2;mzO5K{a!lVL;(rWY4M10Y3W8&#SBha8mO=T8ZZsluWE znFvcdK%l{lVtvew%)WI__}l;@b{5%ewBQw`k@F@NY7F4#?mdQsy(mejizbnS-A_V- z%Bo$pP= zro_<eFd&@f?O^OeC_1li*R>!NvEPl`hOa`+0K6JeH2@nHflI0Z0%b ztW1XP%YNV5jU6}PrJd@qkotseH#JH*C@f|loB)e@n^{j5Mn*BK)|;DKKD}hs#D>lR zRpqL2lSZ7zp_;+HV=vmBOFvEiWbbUeNjInH&Z=-2RlGhi>zmy}*casEbeZ$qC_lK@rN!|t3 zrL1Fv*2y!|8i}E&>xZjX%haC859^<_t9>>4ha7BuN5hmP5~AxG@y@%1q7-7vMllG)$V4&p^)M0g1M921~>C zt8P+@I!`d=Lmx^_2Ty<0Za0Lc{-E&mUZk~;fsukUM1NMilIRU3$??Z9j8=H|!EOWLiMw(MOOM%wM~hND@i^WWVJdwbfWV>`loPT(+tN^w2}aX!5S z)&$%Gk?}*3{M1f|cSqA>JA->&;^U$a1aSkhk*8znp zi6|je^U7N?0E)pWNA!e;elQL`ZyTE>{oAAY!9TEvmCOE|YvTl41Gb4?n!kID{W)Zo zbLyLSUS~NyeK~%vEoYg1M&#ja3!VNd8SLA87IDq(51y@eZnfCL96{|PS(4gmq(k=Z zaVIax(5aTR&|D5LA4hG-cfDe3L}6{)#;=DAJ?l$NjqWRQbuUFkHSs`lI>7KLoNN~f4DR^ADdAUn z<#M!+&4iP_H4pp}kNQ#j)>a6Q0oT~92?&(%n#QJ{u0C! zy(C!t>*Y<~hF?1_+JYC?2v?aG4TIm_D^8kFj1TxL+DfSjG&##BWZ{$gEoUcp_o?l- zz3(@k5{$u~C|M$#YM`&P8F6Q8D=HdKi<4 zzlnFuD834a_iicwwh9RxhbtDt_+pKkdaSWx?t)nK1}tqboCll7x@B&MN*8nO-JO$- zCX493uW&vi64SmiUD~nWg*g-AQR$&4!&r{@b*}=QgLUKKw{>w=KpP6QMG^kFz>6@h z`ii{i9G=l`Vv7&OO;V=8$?^5fEg_P)TTqu2I2l|T%%ms7 zu!EzZr%^K0upt}!k6kS!XQB%<)dRzbA6~S@rBo2RO!a4y1GYPyN`>aDJIo=JTxkA!4Vlc6tXD59FndVAy@z(HF7T7O-i@mg|--ZSiou>?<9~J(-GBU^ZA!uY$HnQEY)bjUDd>?!`o7XO&; zg9HCVdisaF|5I@CvFQIif&GK<4|a0<^{+nRzwGoE&HvJ&{0BWFs=2zD8{3<^$T(Oy zqKiAanb>)e*jYQ+nwx!GK&~Ie#7PeL<@T4D4B!$7mc0BHy{$?OWNMg%VZSqQ)o30u zHHBbl-pf`?*oJ6YVHZNMp5~sc5pwXepI**i1`X^%IGqdc2**U;pq1y)B`>;b><(8F z+@_K_Q&K6KP}7#vsh`W{Jki);Ojagby>Fdet`ujk>BwM}(AA!vWNxl}jJ>~K`M)nI zrP1_|>57YPvm8T=NyK^;(uD3kA#6XSjj6sB0^PgMpb}>F-SErdR)yD^GGg&%l{Q&;j@VYjv>l` zCsnUFc_lt&Hwgst@i|~UA-{FNSe|>HW_@7YmenRZPa6t%K=HL{0PC!sUy|?^HKGrZ z@totlNV}O;PA+;@rUJj0Q1%309jXcHkp_<~QFKHCzzEK{iUGZS_lr8t6L8u|o>0-& z0(f?mgA6rZPHY0k!V)q3)`Uv!)S6*XJ#;xx7Z4kw(kIA%SGG;q31NmgVnE5k37A&}~ z^hD*;Ls>(s$q{|1a2?u_+bNr=&%CRtn|k3M_ZP~Je)vtWTf#f&!Zs3E?g6Z?^O@OM zN3Pgd>k*AbLa^hn+IuP?biQHx{~jODFSISgA_DM^@Y7Gg-GYJ<@p76&IzDskfO5vJ zmm5Blicj%8bedEsXN}KfcG|shjOLTGXvCeT_dP+qx|19*hD>QqqW>1M*w+dZ4(Ns? z`@Y9$^+JG$u}fSN`ow^jBNj|X&Xn29ey#@FXy-G{+uWWoWgKUL@Z0vfdVgOnP9PC9 zt?s!+@Ut1a_8b2jzsH%@efEHIhK;{Qn%HkT6mtVuEU`)eeUM*`9BbL*jg|0<#AgNX z!8-%~gKUgBT?;x(@4Ptg>(|Zn_ub#y{070SoPh$7@hJ!5zFA|gt%6j9HOjfd`siw& zAyj9TvQ&Y{^5rrdJeCGnlv#{&<`xJL0Hp>JM_3qA64M%vmHk5A8PU=Ea#z~^Ef{C6 zDqaMeLnZhj^g#J!<>fkQ#E!oi@9S=g-Ki@>mr zXPoag>+Hy~B6Ur9nY0~c_UKTHHkIsMCo5FR>0wBs5j{Mn;kwa7-htDU45UO1a+vy( z#!)J(rNwh15JNTxizNm^NxulWP0I-+OK^T_Rj8*viQlX-hZo#PgugYe%Cb$${VD_x zFeZg^xx$1NJi@7LKa*~yC@*+(zZ{b?UT7^Rm2GpJ<*}Fm0T&-ROiB&Bp>1o?x$Y4) z*={88>p0^NZ(e%@(onA|gGdYAb#b0tXR*7M&6YjNdLDvnKAsat-XTU+09lJWy%1Cg zo^q#Upmdq^ot`MLvgDQ3cv~xTAnp`XCy8PG#P~Qrs0+kkF4l)7TgThJmKQqr1jYW$ zUMTLLJ(^!8w-|$_fFK$vt(wkIYpokXuzfAA)8dP@A{83=&4W$S+bt%IDmnv2%Fe^F zF2Um~Od*Cd60FOzdviBp09D&r#`k5Id{~5gE#hcJR12p;{`1vN2&r@Z2~tIYTGZzN zmg_SW$>W4LW$%EwjJ+$wRrX|*(BL__Y!>H24tVSFL{*h%7IBDGM$+@2GQy-Sekz3z zhQ`|~mMpVjUBzpCT%hz0)Hcq9w`#$fAmJyU$%6%&PF(u^DQ9`WI#%aT@I>JV6t4bHzfhm zLB};-;;Q#dp))xw45YwOi}<)1)lmy2E`LX^3#ARrnI|}_8{Ds!h8*ZdX!wbI(U%K} zR&r&z(~T!~en$z(jdHwO!yb7+jyF@_*w;cgmMhiLoQq&y+;&ujS}62V6`%7pSYpmh zW&uGC{3|9oJ=>6#o>o%phRcyOm_kP6{mXodqbH& z$dN08p#427dE32eE4iN6f+lVVB-jysojOx%2frLBYtvQAb&kTR)}HL}rbP>iz6bkO zH~w1Yor!B$%{ndwCPkgWYU>teJWCIcWBA?#ro^h|Q#>z=6_mDcqif25l>Rxs@nMKu zK1J4Cz9EsoXAYXeIY>=H2rGCtx9-{`dD}{lGX-PkX}Qj}8<8!cgc<&+lc5%0o+%8< ztxDNWORNk0!bT`5VcwxZJcBf#Y@jUkJ>Q z7|btWouF(ZCn;S0=@$c0h^{GjC4PtU?r8}7p;qFjcvQ^Z4T$An-(_=DI#kT^%SMkFK3C=oc5|1bRdDA`UTv|q3|r?QqDlr)?g~y@h;-`C z)$SKjoTrq5%#4-wEMD2D#=DA%XEX$x$r&u;d|2yeAB>eOsi`NDyjDE$q38Cy)o=LI z&O|de(cpIpCRx;oF$#Du5PXYiW{<0vR#bbuZ+czc?e;&G)wYT8S2yz-4A?o%yRMBJ zF86aLz@4X+shO7Iu`*DYbP_@EFSs8&*E%gz&+R%I4p`sfm^kMZ+QZS$dgYbclOH&B z;+%}!&DhKwUwP@7n>akUU%+o7A~{kLlwv0BvN;wc3R~q>dXlNJ2embuY045*m}i-D zzMs-3tegj4+6D3BL4+}>;CXq>$jeeznaj@D5LMfAU0!ZuG0vI`%H7=!u_j~* zyfE3nAV2#@$A;_7)GpaH=c22Lt|fG^`f$xLdZ-Hlgw{Ky>&E#l34_Ey^+1ZjE+RG{ z1#4@)E9JnWk>se<#`}=O&i7=Ywy56hu{qX$5G?&mivkix`l*zebzglHIe$jE6YxGE z@SBXq`H$U@O$3imv%6ochx_B@EqhF6^Gxd#yT@}t%M%z8hnm zm+`h5G@<)(fF65B3cIgC;YG-aEU4wQTL;#RhBU7IO!iok@V02$D)A*_c6^cM8VPSt zbXP4ttw~>-) zHE>;mBT{2+kr<}5PueB@wE~4ZP74&ZLY9gwn!-l-vtNS^C=||!Y|7)qrcD8U`YOUG z)FYZ3=F;f9{k-a$teN_x9Q_vy*_mf^o`)12mGLi=q@-D0tSEeM^XR`d7Thh)rw zLIO9zvI@W&hbcnCq?bD$k}s|nRSD|1Mov@i#R$vA8N#3__Rv(NO_*;fO3J&4t%WDV zJ*oeA@T6R|IgxQ9hEHJj>yX4UkWzj3j6;CvGm~UZ>v7Q9BTe0^3oGXEi_VbRYZu;I zO+iM>r)Q@gc4VP0Vgn?!YDfwt;CH_@AF~854|69ih0}K?2t%AQ8H6IRHgAj*bkn&Y z*;P@&Ci5K^6}G|$=wheiUM)lJ86-kc_NfsNPhijQq06spwZmmXgauORrz@{{Dm2<9SUJIer|}vz~$wp&AKDm5kk{ygOcU~4#m_Mut)vd7J=1e z(MJC(Uv&VdJf_%px%t>>c=I&~tV2gJ$8z&dqwd`4_agb_mn=OUgE){zfDUa)nKc<#` zW1~lWalyr|vA8wf)!^S_lV>XOUqoF%;eLiB>;(K=HREf3J!~`3^ZwIzxo-aW_^{R8 z3usig7|o|6@6Nc|Em!B$7}aH*iQHEd=CA!s2|IdEo6(whNx#$T(G zq{Xy=L-#s*DQdQ!G*|Qbd~vN=2gRX`C>QC<`Pg1h>CFOE9EoG3)Qj@@wrBqrCm*MX z#HdL5E*4=pB`zNJ4|dz9WKT&bT;90OMKT)C1)Pp$peO>cc>ClQO^x?LRS@UAWaU&QBIH?5UIBjSd8Wy1iS(NiH%k1#WGRsG&*HRspTK2Dh(C7N#OP zQq*~?`MqczUDMH ztX8pdzO4a{3_X>SlHh%8DbEOCR3tbslM4;RCaE~}Al24t?XH(5=>-$&ma^8*Jxy7DUi#jO-RaQIEmX(*MvFouCI+(@P@uLX~ zR07#>WTh0~up6+@a~zPR3JyI@c?hbbh8pFi zhQ>?i-k(%8WYU5o2Fj#;OODBtj~Ol6uL^PW^vVc3oa3fVxr^Jo5PD7YQ{(nRe>ImW z>!mxaM@&IMSOS}=apx${TvCr#)uzS0Lo{XbjnBX@AH32Eht7?mTxRxAZ+v#b+(Bii zNw6tpfQRU#7|wW~&oU-6h+MiY!JR)NX-XT%LY?5+tEq?Uv_{FIL_I-agDn_(Y@*nr zqS|^WRmReU4SU_54TRdo<;+t!GNerBti*$q4$2=^l(3HAVd2 z7os8gy(n5mvLAygmOe5sIK)&RhiUbxI>k06MH%DA4$7FMVjF`DtF zSx;?7Bse0wAB+A}r26T`w^J9O8J{{A?byrlM+`jV2y2J8X#K3U!W|1~wooK5GTCC- zG?iQ*#6*8zj}(e}klHz;&z&r3rJnXRnr`gzsPC(GA0cd!F=qiq1qm+dfj3w#>G_vy z^xq!7J$OJDtt?`CcBE`#V?4_x@tU(Ip;PQQ<`bN9Aqrg zb{w{%pP`kMa03^h8ml~F<|mgQBlzg_X@5arpxhNOqb@hUgfh^={xuW7_tOFmhqt?% z(U3l5ON{RNeT6B5czecJSV~m8DhSugF2#1Bp?@Sym9RIcq;Crk<7t@!;qk%y0yC(> zET~MjE?N{zBRfoGk}WQD9%BH!PHO?-b(gk^Vzc;RW&1!E_c%ou2=}Q0+L;S)pP9}` zZLlFh{7{FGzS~y8g}{vv3el_~SLzEZ3$aRwz!Xa~B8Wv@T{I$&A|lhJzc;blf04x( zrHO`StWX$=j6%K9+(Ly-f|9KwcLan~t|h3EW;iB=8*GUi{Kpkhle(uU8Qe8Y-AM%_ zeM0dPTry=6PK6}x?z1dB@fUcnNX*6Y@MLtE04T__SsK-EF5~FBDC#Od^4O}2+Wa&n zZnWaXVo6ij`hPl}4MWU{+eg!L&a?GP#XE-me84}W1Xj_ZJ<#=Nhw6HxfAXMdF4AAj zDqvIEgda2}p)|LQ->vb|M|kV+;5j#Y#FF*u-q`-#9?I~;-q{WEDEy2hY0|t<$s}T2 zFC8=S8%>mtE{^coiYH&rT=~$a3E0F&+OA_n&r-N)-bsw{zKWF<{VzZ8GnKkQN;3(B zZm+z&644HYR!X-mpVEWAAc#K0enR*b-@5A zP5o|z$k$pfS68q?f(CE>Nh_Or!s}~sm0E+XVRRv7rpUlV+5o|(m1re!675`l^2W1&m>9G$ zU^~%yiKMSWs8A8c7$(FZ7roOgAF^deX6E*3EzK34mJM>Oh33Gljs+|_^ol$^tDRa! zT|BHPR|@xs-Tc5TmE{T4%?X*Ic$qJ}s#Auvui`W*dmqW&4qbx7=eYx#bd+XeS&LyG zM^M;ova+8F_&0sMG$P}Rc3#F79r?c8M0$;d0vbx;qge zK-w~s@&;CJXbmGWT8zi{n4w^&qI$Ijp-YmOl*j<4;+-CK3$iHeYM$)9@Sot3H{%Jb zPvc57ZSpa{A@c$vsRi|ASn(j;eoG#W1&Hm2Hl-kSepbJ~1vA^5#7m+ReB13)5uWMa z)26YGf$JWd=tdLw=DP;QJ2x>5p`{xaSJ<$wD+A#2+srmIAFjaCo#tXKo>kA3CF1)E zhMx@>tZzA=Sh8e`?G;T#hTPsXJYZ^l=w$jTp`CAZMfG_bRxyeojI`)$XtO43^DSt* zu{bP^iz>Sb1Wx@Tan@g^zYHSG5D%LgsnPGgI&-PK5|Z;uV@dNwtD{koh@K9YKBFDP;03#`PSp0pdE+BYLx~PP z;&tDi{jeo~lYsDaAJh%y??s_P|HBj4mCQ(;k~YZvxFL#MFdq;@*#STBL!nHJEA@U` zZ6E)p3=64#B%YT18O>NbW!!$$lEJcrhxg)1E!6dDZj|*Lkh1&gUsMqILS8@ABNoW~ zjk#r@*Wo#&0+is}gN<-+HjoSPJ&g4v8#}poDr927iUa0gXh4pVAxsl3H2aoPr_v9}P4%ssK&7N;>7LE@ZD=XbYB>D!^M-CH}g z`IGJRO%K?u4lfC6e8i3taH+|Q&Rrp$j?B}AO;sN&8Rf_15^tKxZDA2K{|*`r*x9z2 zxwRQ^qFOWMbfl{EoBGY<1f_p)T=sZDKe?H9PKACWVvya0hvw`~l!xXu?aG0+@i1Kh z7$H#wSf<lb50|vRCRJ)M|b#j|S z&aKVNH_c^6^%n;acdQwPA{!05nE6E3j||TFKV(J<8sP8}@G-#oFw00OGx+W1y-R4d0 zp@F-~2&0v)z^0Ib;ui2K?7&_mdND(Z(oTOABb7=${8;1|_oUgN_s=T(_H`>wC>RMVK@y zF3x6^exwu&Wm&$^8L{gz0Em|YG#O^kvttyP>dRy8n?m3(3|ka9zvLE}PD@Ej?h|+= z7o4s3SlY`yVLz}PF>f!ko`s}WveB_bP>9Mt2S*8>C==Cg^rtPz`i{tcS5L3z>>08>q&uR~Qkb$UI`U%j5JwwmMa z1OFpVzV=|F$}Lo#PKRhLS-4)T+rLH*dN(q1fx>RZ`EP_$LwjMGS@X1Xw- zBZo;o2NCUtzh4P$)^EWr(+ure5reG7=jwZ65g6%I+appxu<4oEA+_MSmg=VXcVU-zn#Ll zLf&$?6|m~0&ZO>JQ%YlT-D%zFY#GvMDhBx$U+2*o?i+TQ;`g=5)z><(CAI~+xqyf_ z$maEcr0Ra5ey-8qFI58J2xmSGAeiqsM2u?`CUVn%^Fx8 zPdOgVEK892^_LQQ0~^;xx*i72D1rd=69_mhO-NRv^;-=yYN>)9ogC9Hb8&81KY!Ed zM}D$(3xO4(s4LjtV8p{ZcG+Pm$DCFh1VG;+y^p_&r>`awtc@S=dyaj6`3yPUIuT1n zTcHpnPaT^ainnqd5g)(*T0aH!9?rizhr^J-s>#twen+#RBrCe*d6qsEc{MCJ zv6g!-D>bh~*Y#9hl9eWp5Pg6&i|rc$;bh;mPByZV)=pIci6Q?Tq%zRJ5-hp+OqxB6 z!VzKJGD;b&o6~_5r5hFCjS}~fi3OYB&B1n1SB=LM05+={`j6-< z^A#0)P0SMSdU2HNk%`5bL!|^U0x3+oWR1~!iD;td@yd|?-3WZdrX+B84uYz8O7PyEQiNCPScT8BS-lsK5uX4+n^9c*$P?3rNHVwbzFh>gCgmE zf$qv<0CBf=F7hd}#&N%YCk>wy3L2rIX4EaQn-_`1Zxe$|*jr67d{YGZdnx_gt2{`XB6rN;v8wYwnz$VXgB6pj7C|z(-@`f^` z`nc9Df`A!W=IajnNP|q}W;%Cqu&%0S#uv!KMNk7>Q>7XFK|uRNFIj|uauo(jK!&!P z%v!~m6Zvu}R~1Ip#NmeWIgD|=t0PG(;pT`glE{T%eYmHmpIqa<0wpzS7K*UV>qx$@ z_UYGQ56t{NgNdVhG+di?RBw$#3ZwaMpUcIk>1Y*on(*F;3W2VfsmaHHjwK#{(wxXf zuB2?QUJdCxAVD`O&W`HaXW5Z=Wy^=D^t73tY&S_M_Q8~gZ-E^~l4rD?2pK~bH4yTt z{o;Ib)#GVay)?Y<5lMID4+F`9-X8Yg-y^<42S=FoS`3|jK{kt|Ly~7~!o5el72lcV zaV~=43?0=(SzAb5cK<$=HZ9NB)W&J52PNm2T2#v5(5$xsZ9EaMDM}f^(y7qaS!}3O zVUd-UxJ+fzB+G5mpke_Yr6K@kf6rnoHbHzNLz&Z8baxcF-Yt@^l0#CUH=GzHq;({A zjgy@1)9DmSQZ^SGxHuFHGy0e)MEV-uPRAN$gR_&Ch~uXJMV;P4b=WKU)f>tC7M=X3 z$jy$;pl~hSh_FleF0WHp8G%nEf@piU4k>g3xTMIUX46KPUTGWJJHi-ZJm9|oSYrAi diff --git a/roofit/roofitcore/test/TestStatistics/testPlot.cxx b/roofit/roofitcore/test/TestStatistics/testPlot.cxx index fbedc0968468b..8a3ee13948659 100644 --- a/roofit/roofitcore/test/TestStatistics/testPlot.cxx +++ b/roofit/roofitcore/test/TestStatistics/testPlot.cxx @@ -11,83 +11,99 @@ */ #include +#include #include -#include -#include -#include +#include #include #include #include #include -#include #include -#include - #include +#include #include using namespace RooFit; -class TestRooRealLPlot : public RooUnitTest { -public: - TestRooRealLPlot(TFile &refFile, bool writeRef, int verbose) - : RooUnitTest("Plotting and minimization with RooFit::TestStatistics", &refFile, writeRef, verbose){}; - bool testCode() override - { - - // C r e a t e m o d e l a n d d a t a - // --------------------------------------- - // Constructing a workspace with pdf and dataset - RooWorkspace w("w"); - w.factory("expr::Nexp('mu*S+B',mu[1,-1,10],S[10],B[20])"); - w.factory("Poisson::model(Nobs[0,100],Nexp)"); - w.var("Nobs")->setBins(4); - RooDataSet d("d", "d", *w.var("Nobs")); - w.var("Nobs")->setVal(25); - d.add(*w.var("Nobs")); - - // P e r f o r m a p a r a l l e l l i k e l i h o o d m i n i m i z a t i o n - // -------------------------------------------------------------------------------- - - // Creating a RooAbsL likelihood - std::unique_ptr likelihood{w.pdf("model")->createNLL(d, ModularL(true))}; - - // Creating a minimizer and explicitly setting type of parallelization - std::size_t nWorkers = 1; - RooMinimizer::Config cfg; - cfg.parallelize = nWorkers; - cfg.enableParallelDescent = false; - cfg.enableParallelGradient = true; - RooMinimizer m(*likelihood, cfg); - - // Minimize - m.migrad(); - - // C o n v e r t t o R o o R e a l L a n d p l o t - // --------------------------------------------------- - RooPlot *xframe = w.var("mu")->frame(-1, 10); - likelihood->plotOn(xframe, RooFit::Precision(1)); - - // --- Post processing for RooUnitTest --- - regPlot(xframe, "TestRooRealLPlot_plot"); - - return true; - } -}; - +/// Plotting and minimization with RooFit::TestStatistics: minimize a +/// RooFit::TestStatistics likelihood with the parallel gradient and plot it as +/// a function of the parameter. The results are validated against analytic +/// expectations and against direct evaluations of the likelihood, instead of +/// the RooUnitTest reference file that was used before. TEST(TestStatisticsPlot, RooRealL) { - // Run the RooUnitTest and assert that it succeeds with gtest - - RooUnitTest::setMemDir(gDirectory); - - gErrorIgnoreLevel = kWarning; - - TFile fref("TestStatistics_ref.root"); + RooHelpers::LocalChangeMsgLevel changeMsgLvl{RooFit::WARNING}; + + // C r e a t e m o d e l a n d d a t a + // --------------------------------------- + // Constructing a workspace with pdf and dataset + RooWorkspace w("w"); + w.factory("expr::Nexp('mu*S+B',mu[1,-1,10],S[10],B[20])"); + w.factory("Poisson::model(Nobs[0,100],Nexp)"); + w.var("Nobs")->setBins(4); + RooDataSet d("d", "d", *w.var("Nobs")); + w.var("Nobs")->setVal(25); + d.add(*w.var("Nobs")); + + RooRealVar &mu = *w.var("mu"); + + // P e r f o r m a p a r a l l e l l i k e l i h o o d m i n i m i z a t i o n + // -------------------------------------------------------------------------------- + + // Creating a RooAbsL likelihood + std::unique_ptr likelihood{w.pdf("model")->createNLL(d, ModularL(true))}; + + // Creating a minimizer and explicitly setting type of parallelization + std::size_t nWorkers = 1; + RooMinimizer::Config cfg; + cfg.parallelize = nWorkers; + cfg.enableParallelDescent = false; + cfg.enableParallelGradient = true; + RooMinimizer m(*likelihood, cfg); + + // Minimize + m.setPrintLevel(-1); + m.migrad(); + + // The analytic maximum likelihood estimate is at Nexp = Nobs, i.e. + // mu = (Nobs - B) / S. The tolerance is at the scale of the Minuit + // convergence criterion, given that sigma(mu) = sqrt(Nobs) / S = 0.5. + EXPECT_NEAR(mu.getVal(), 0.5, 0.05 * 0.5); + + // C o n v e r t t o R o o R e a l L a n d p l o t + // --------------------------------------------------- + std::unique_ptr xframe{mu.frame(-1, 10)}; + likelihood->plotOn(xframe.get(), RooFit::Precision(1)); + RooCurve *curve = xframe->getCurve(); + ASSERT_NE(curve, nullptr); + ASSERT_GT(curve->GetN(), 1); + + // Every point of the plotted curve must match a direct evaluation of the + // likelihood at that parameter value + for (int i = 0; i < curve->GetN(); ++i) { + const double muVal = curve->GetPointX(i); + if (muVal < mu.getMin() || muVal > mu.getMax()) + continue; + mu.setVal(muVal); + const double directVal = likelihood->getVal(); + EXPECT_NEAR(curve->GetPointY(i), directVal, 1e-6 * std::max(1.0, std::abs(directVal))) + << "curve point " << i << " at mu = " << muVal; + } - TestRooRealLPlot plotTest{fref, false, 0}; - bool result = plotTest.runTest(); - ASSERT_TRUE(result); + // The likelihood must also match the analytically known Poisson -log L, + // which is defined up to a mu-independent constant. The comparison is + // restricted to moderate Nexp values, where the truncation of the Poisson + // normalization to the observable range is negligible. + auto analyticNll = [](double muVal) { + const double nexp = 10 * muVal + 20; + return nexp - 25 * std::log(nexp); + }; + mu.setVal(0.5); + const double nllOffset = likelihood->getVal() - analyticNll(0.5); + for (double muVal : {0.0, 1.0, 2.0, 3.0}) { + mu.setVal(muVal); + EXPECT_NEAR(likelihood->getVal(), analyticNll(muVal) + nllOffset, 1e-4) << "at mu = " << muVal; + } }