From ca77f818afb9dcdd4ecd2ce951f6388094760b74 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 6 Jul 2023 16:06:33 +0200 Subject: [PATCH 1/5] [RF] Add new version of the BarlowBeeston tutorial in HistFactory style --- .../roofit/rf709b_BarlowBeestonHistFactory.C | 252 ++++++++++++++++++ 1 file changed, 252 insertions(+) create mode 100644 tutorials/roofit/rf709b_BarlowBeestonHistFactory.C diff --git a/tutorials/roofit/rf709b_BarlowBeestonHistFactory.C b/tutorials/roofit/rf709b_BarlowBeestonHistFactory.C new file mode 100644 index 0000000000000..d1840ec261cdf --- /dev/null +++ b/tutorials/roofit/rf709b_BarlowBeestonHistFactory.C @@ -0,0 +1,252 @@ +/// \file +/// \ingroup tutorial_roofit +/// \notebook -js +/// Implementing the Barlow-Beeston method for taking into account the statistical +/// uncertainty of a Monte-Carlo fit template. +/// +/// \macro_image +/// \macro_code +/// \macro_output +/// +/// New implementation of the original rf709_BarlowBeeston tutorial, using the +/// same kind of RooFit models that are built by HistFactory. Therefore, this +/// tutorial is a good starting point if you want to understand how HistFactory +/// works. +// +/// \date June 2023 +/// \author Jonas Rembser, CERN + +#include "RooRealVar.h" +#include "RooGaussian.h" +#include "RooUniform.h" +#include "RooDataSet.h" +#include "RooDataHist.h" +#include "RooHistFunc.h" +#include "RooRealSumPdf.h" +#include "RooParamHistFunc.h" +#include "RooHistConstraint.h" +#include "RooProdPdf.h" +#include "RooPlot.h" + +#include "TCanvas.h" +#include "TPaveText.h" + +#include +#include + +RooArgSet createConstraints(RooArgList const &gammas, RooArgList const &samples, RooArgSet &globalObservables) +{ + RooArgSet constraints; + for (std::size_t i = 0; i < gammas.size(); ++i) { + RooAbsReal *gamma = static_cast(&gammas[i]); + + auto nomName = "nom_" + std::string(gamma->GetName()); + auto tauName = std::string(gamma->GetName()) + "_tau"; + auto meanName = std::string(gamma->GetName()) + "_mean"; + auto poissonName = std::string(gamma->GetName()) + "_constraint"; + + double nomVal = 0.0; + for (auto *func : static_range_cast(samples)) { + nomVal += func->dataHist().weight(i); + } + + auto nominal = std::make_unique(nomName.c_str(), nomName.c_str(), nomVal, 0.0, RooNumber::infinity()); + nominal->setConstant(true); + globalObservables.add(*nominal); + + auto tau = std::make_unique(tauName.c_str(), tauName.c_str(), nominal->getVal()); + + auto mean = std::make_unique(meanName.c_str(), meanName.c_str(), *gamma, *tau); + mean->addOwnedComponents(std::move(tau)); + + auto poisson = std::make_unique(poissonName.c_str(), poissonName.c_str(), *nominal, *mean); + poisson->addOwnedComponents(std::move(nominal)); + poisson->addOwnedComponents(std::move(mean)); + + constraints.addOwned(std::move(poisson)); + } + return constraints; +} + +void rf709b_BarlowBeestonHistFactory() +{ + using namespace RooFit; + + RooMsgService::instance().getStream(1).removeTopic(NumIntegration); + + // First, construct a likelihood model with a Gaussian signal on top of a uniform background + RooRealVar x("x", "x", -20, 20); + x.setBins(25); + + RooRealVar meanG("meanG", "meanG", 1, -10, 10); + RooRealVar sigG("sigG", "sigG", 1.5, -10, 10); + RooGaussian g("g", "Gauss", x, meanG, sigG); + RooUniform u("u", "Uniform", x); + + // Generate the data to be fitted + std::unique_ptr sigData(g.generate(x, 50)); + std::unique_ptr bkgData(u.generate(x, 1000)); + + RooDataSet sumData("sumData", "Gauss + Uniform", x); + sumData.append(*sigData); + sumData.append(*bkgData); + + // Make histogram templates for signal and background. + // Let's take a signal distribution with low statistics and a more accurate + // background distribution. + // Normally, these come from Monte Carlo simulations, but we will just generate them. + std::unique_ptr dh_sig(g.generateBinned(x, 50)); + std::unique_ptr dh_bkg(u.generateBinned(x, 10000)); + + // ***** Case 0 - 'Rigid templates' ***** + + // Construct histogram shapes for signal and background + RooHistFunc p_h_sig("p_h_sig", "p_h_sig", x, *dh_sig); + RooHistFunc p_h_bkg("p_h_bkg", "p_h_bkg", x, *dh_bkg); + + // Construct scale factors for adding the two distributions + RooRealVar Asig0("Asig", "Asig", 1, 0.01, 5000); + RooRealVar Abkg0("Abkg", "Abkg", 1, 0.01, 5000); + + // Construct the sum model + RooRealSumPdf model0("model0", "model0", {p_h_sig, p_h_bkg}, {Asig0, Abkg0}, true); + + // ***** Case 1 - 'Barlow Beeston' ***** + + // Construct parameterized histogram shapes for signal and background + RooArgList gammas_sig1{ParamHistFunc::createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10)}; + RooArgList gammas_bkg1{ParamHistFunc::createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10)}; + + ParamHistFunc p_phf_sig1{"p_phf_sig1", "p_phf_sig1", x, gammas_sig1}; + ParamHistFunc p_phf_bkg1{"p_phf_bkg1", "p_phf_bkg1", x, gammas_bkg1}; + + RooProduct p_p_sig1{"p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1}; + RooProduct p_p_bkg1{"p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1}; + + RooRealVar Asig1("Asig", "Asig", 1, 0.01, 5000); + RooRealVar Abkg1("Abkg", "Abkg", 1, 0.01, 5000); + + // Construct the sum of these + RooRealSumPdf model1_tmp("sp_ph", "sp_ph", {p_p_sig1, p_p_bkg1}, {Asig1, Abkg1}, true); + + // Construct the subsidiary poisson measurements constraining the histogram parameters + // These ensure that the bin contents of the histograms are only allowed to vary within + // the statistical uncertainty of the Monte Carlo. + RooArgSet globs1; + RooArgSet constraints_sig1{createConstraints(gammas_sig1, p_h_sig, globs1)}; + RooArgSet constraints_bkg1{createConstraints(gammas_bkg1, p_h_bkg, globs1)}; + + // Construct the joint model with template PDFs and constraints + RooProdPdf model1("model1", "model1", {model1_tmp, constraints_sig1, constraints_bkg1}); + + // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** + + // Construct the histogram shapes, using the same parameters for signal and background + // This requires passing the first histogram to the second, so that their common parameters + // can be re-used. + // The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`. + // This allows bin 0 to fluctuate up and down. + // Then, the SAME parameters are connected to the background histogram, so the bins flucutate + // synchronously. This reduces the number of parameters. + RooArgList gammas2{ParamHistFunc::createParamSet("p_ph_gamma", x.numBins(), 0, 10)}; + + ParamHistFunc p_phf2{"p_phf2", "p_phf2", x, gammas2}; + + RooProduct p_p_sig2{"p_p_sig2", "p_p_sig2", p_h_sig, p_phf2}; + RooProduct p_p_bkg2{"p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2}; + + RooRealVar Asig2("Asig", "Asig", 1, 0.01, 5000); + RooRealVar Abkg2("Abkg", "Abkg", 1, 0.01, 5000); + + // As before, construct the sum of signal2 and background2 + RooRealSumPdf model2_tmp("sp_ph", "sp_ph", {p_p_sig2, p_p_bkg2}, {Asig2, Abkg2}, true); + + // Construct the subsidiary poisson measurements constraining the statistical fluctuations + RooArgSet globs2; + RooArgSet constraints2{createConstraints(gammas2, {p_h_sig, p_h_bkg}, globs2)}; + + // Construct the joint model + RooProdPdf model2("model2", "model2", {model2_tmp, constraints2}); + + // ************ Fit all models to data and plot ********************* + + std::unique_ptr result0{model0.fitTo(sumData, PrintLevel(0), Save())}; + std::unique_ptr result1{model1.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs1))}; + std::unique_ptr result2{model2.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs2))}; + + TCanvas *can = new TCanvas("can", "", 1500, 600); + can->Divide(3, 1); + + TPaveText pt(-19.5, 1, -2, 25); + pt.SetFillStyle(0); + pt.SetBorderSize(0); + + can->cd(1); + auto frame = x.frame(Title("No template uncertainties")); + // Plot data to enable automatic determination of model0 normalisation: + sumData.plotOn(frame); + model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0)); + // Plot data again to show it on top of model0 error bands: + sumData.plotOn(frame); + // Plot model components + model0.plotOn(frame, LineColor(kBlue)); + model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure)); + model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed)); + model0.paramOn(frame); + + sigData->plotOn(frame, MarkerColor(kBlue)); + frame->Draw(); + + for (auto text : {"No template uncertainties", "are taken into account.", "This leads to low errors", + "for the parameters A, since", "the only source of errors", "are the data statistics."}) { + pt.AddText(text); + } + pt.DrawClone(); + + can->cd(2); + frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately")); + sumData.plotOn(frame); + model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1)); + // Plot data again to show it on top of error bands: + sumData.plotOn(frame); + model1.plotOn(frame, LineColor(kBlue)); + model1.plotOn(frame, Components(p_p_sig1), LineColor(kAzure)); + model1.plotOn(frame, Components(p_p_bkg1), LineColor(kRed)); + model1.paramOn(frame, Parameters({Asig1, Abkg1})); + + sigData->plotOn(frame, MarkerColor(kBlue)); + frame->Draw(); + + pt.Clear(); + for (auto text : {"With gamma parameters, the", "signal & background templates", "can adapt to the data.", + "Note how the blue signal", "template changes its shape.", "This leads to higher errors", + "of the scale parameters A."}) { + pt.AddText(text); + } + pt.DrawClone(); + + can->cd(3); + frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)")); + sumData.plotOn(frame); + model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2)); + // Plot data again to show it on top of model0 error bands: + sumData.plotOn(frame); + model2.plotOn(frame, LineColor(kBlue)); + model2.plotOn(frame, Components(p_p_sig2), LineColor(kAzure)); + model2.plotOn(frame, Components(p_p_bkg2), LineColor(kRed)); + model2.paramOn(frame, Parameters({Asig2, Abkg2})); + + sigData->plotOn(frame, MarkerColor(kBlue)); + frame->Draw(); + + pt.Clear(); + for (auto text : {"When signal and background", "template share one gamma para-", "meter per bin, they adapt less.", + "The errors of the A parameters", "also shrink slightly."}) { + pt.AddText(text); + } + pt.DrawClone(); + + std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl; + std::cout << "Asig [BB ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl; + std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl; +} From d0ceebb50a5a171ca01d100f2892c31c73bdc3f3 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 6 Nov 2026 10:46:16 +0000 Subject: [PATCH 2/5] [RF] Fix empty-bin gammas in rf709b and move it next to rf709 Gammas of empty template bins were floating and unconstrained, i.e. flat directions of the likelihood, so HESSE failed for the Barlow-Beeston case. Fix them like HistFactory does. Also move the tutorial to tutorials/roofit/roofit with the matching \ingroup and index entry, add the includes needed for ACLiC, and update stale comments. Assisted-by: claude-code:claude-opus-5-5 --- tutorials/roofit/roofit/index.md | 1 + .../rf709b_BarlowBeestonHistFactory.C | 33 ++++++++++++------- 2 files changed, 23 insertions(+), 11 deletions(-) rename tutorials/roofit/{ => roofit}/rf709b_BarlowBeestonHistFactory.C (91%) diff --git a/tutorials/roofit/roofit/index.md b/tutorials/roofit/roofit/index.md index 347248a7f765f..09e7ee550e84e 100644 --- a/tutorials/roofit/roofit/index.md +++ b/tutorials/roofit/roofit/index.md @@ -180,6 +180,7 @@ Explore the tutorials below to discover the main features of RooFit. A more inde | | rf617_simulation_based_inference_multidimensional.py | Use Simulation Based Inference (SBI) in multiple dimensions in RooFit.| | | rf618_mixture_models.py | Use of mixture models in RooFit.| | rf709_BarlowBeeston.C| rf709_BarlowBeeston.py | Implementing the Barlow-Beeston method for taking into account the statistical uncertainty of a Monte-Carlo fit template.| +| rf709b_BarlowBeestonHistFactory.C| | Implementing the Barlow-Beeston method with the same kind of RooFit models that HistFactory builds.| | rf710_roopoly.C| rf710_roopoly.py | Taylor expansion of RooFit functions using the taylorExpand function with RooPolyFunc.| | rf711_lagrangianmorph.C| rf711_lagrangianmorph.py | Morphing effective field theory distributions with RooLagrangianMorphFunc. A morphing function as a function of one coefficient is setup and can be used to obtain the distribution for any value of the coefficient.| | rf712_lagrangianmorphfit.C| rf712_lagrangianmorphfit.py | Performing a simple fit with RooLagrangianMorphFunc. A morphing function is setup as a function of three variables and a fit is performed on a pseudo-dataset.| diff --git a/tutorials/roofit/rf709b_BarlowBeestonHistFactory.C b/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.C similarity index 91% rename from tutorials/roofit/rf709b_BarlowBeestonHistFactory.C rename to tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.C index d1840ec261cdf..0c58cf5bb1eab 100644 --- a/tutorials/roofit/rf709b_BarlowBeestonHistFactory.C +++ b/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.C @@ -1,5 +1,5 @@ /// \file -/// \ingroup tutorial_roofit +/// \ingroup tutorial_roofit_main /// \notebook -js /// Implementing the Barlow-Beeston method for taking into account the statistical /// uncertainty of a Monte-Carlo fit template. @@ -23,14 +23,19 @@ #include "RooDataHist.h" #include "RooHistFunc.h" #include "RooRealSumPdf.h" -#include "RooParamHistFunc.h" -#include "RooHistConstraint.h" +#include "RooProduct.h" +#include "RooConstVar.h" +#include "RooPoisson.h" #include "RooProdPdf.h" +#include "RooFitResult.h" #include "RooPlot.h" +#include "RooMsgService.h" +#include "RooStats/HistFactory/ParamHistFunc.h" #include "TCanvas.h" #include "TPaveText.h" +#include #include #include @@ -38,7 +43,7 @@ RooArgSet createConstraints(RooArgList const &gammas, RooArgList const &samples, { RooArgSet constraints; for (std::size_t i = 0; i < gammas.size(); ++i) { - RooAbsReal *gamma = static_cast(&gammas[i]); + auto *gamma = static_cast(&gammas[i]); auto nomName = "nom_" + std::string(gamma->GetName()); auto tauName = std::string(gamma->GetName()) + "_tau"; @@ -50,6 +55,15 @@ RooArgSet createConstraints(RooArgList const &gammas, RooArgList const &samples, nomVal += func->dataHist().weight(i); } + // Like in HistFactory, the gamma of an empty bin is fixed and not + // constrained. Otherwise, it would be a flat direction of the likelihood. + // The other gammas get the relative Poisson uncertainty as initial error. + if (nomVal == 0.0) { + gamma->setConstant(true); + continue; + } + gamma->setError(1.0 / std::sqrt(nomVal)); + auto nominal = std::make_unique(nomName.c_str(), nomName.c_str(), nomVal, 0.0, RooNumber::infinity()); nominal->setConstant(true); globalObservables.add(*nominal); @@ -141,13 +155,10 @@ void rf709b_BarlowBeestonHistFactory() // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** - // Construct the histogram shapes, using the same parameters for signal and background - // This requires passing the first histogram to the second, so that their common parameters - // can be re-used. - // The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`. - // This allows bin 0 to fluctuate up and down. - // Then, the SAME parameters are connected to the background histogram, so the bins flucutate - // synchronously. This reduces the number of parameters. + // Construct the histogram shapes, using the same parameters for signal and background. + // There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to + // fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background + // histogram, so the bins fluctuate synchronously. This reduces the number of parameters. RooArgList gammas2{ParamHistFunc::createParamSet("p_ph_gamma", x.numBins(), 0, 10)}; ParamHistFunc p_phf2{"p_phf2", "p_phf2", x, gammas2}; From a4bb6d9842b147a8e65eb1e0d2e66342e98cacbf Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 6 Nov 2026 11:01:42 +0000 Subject: [PATCH 3/5] [RF] Add Python version of rf709b_BarlowBeestonHistFactory Assisted-by: claude-code:claude-opus-5-5 --- tutorials/roofit/roofit/index.md | 2 +- .../roofit/rf709b_BarlowBeestonHistFactory.py | 251 ++++++++++++++++++ 2 files changed, 252 insertions(+), 1 deletion(-) create mode 100644 tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py diff --git a/tutorials/roofit/roofit/index.md b/tutorials/roofit/roofit/index.md index 09e7ee550e84e..f4ea1335e6558 100644 --- a/tutorials/roofit/roofit/index.md +++ b/tutorials/roofit/roofit/index.md @@ -180,7 +180,7 @@ Explore the tutorials below to discover the main features of RooFit. A more inde | | rf617_simulation_based_inference_multidimensional.py | Use Simulation Based Inference (SBI) in multiple dimensions in RooFit.| | | rf618_mixture_models.py | Use of mixture models in RooFit.| | rf709_BarlowBeeston.C| rf709_BarlowBeeston.py | Implementing the Barlow-Beeston method for taking into account the statistical uncertainty of a Monte-Carlo fit template.| -| rf709b_BarlowBeestonHistFactory.C| | Implementing the Barlow-Beeston method with the same kind of RooFit models that HistFactory builds.| +| rf709b_BarlowBeestonHistFactory.C| rf709b_BarlowBeestonHistFactory.py | Implementing the Barlow-Beeston method with the same kind of RooFit models that HistFactory builds.| | rf710_roopoly.C| rf710_roopoly.py | Taylor expansion of RooFit functions using the taylorExpand function with RooPolyFunc.| | rf711_lagrangianmorph.C| rf711_lagrangianmorph.py | Morphing effective field theory distributions with RooLagrangianMorphFunc. A morphing function as a function of one coefficient is setup and can be used to obtain the distribution for any value of the coefficient.| | rf712_lagrangianmorphfit.C| rf712_lagrangianmorphfit.py | Performing a simple fit with RooLagrangianMorphFunc. A morphing function is setup as a function of three variables and a fit is performed on a pseudo-dataset.| diff --git a/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py b/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py new file mode 100644 index 0000000000000..a06d6876944cc --- /dev/null +++ b/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py @@ -0,0 +1,251 @@ +## \file +## \ingroup tutorial_roofit_main +## \notebook -js +## Implementing the Barlow-Beeston method for taking into account the statistical +## uncertainty of a Monte-Carlo fit template. +## +## \macro_image +## \macro_code +## \macro_output +## +## New implementation of the original rf709_BarlowBeeston tutorial, using the +## same kind of RooFit models that are built by HistFactory. Therefore, this +## tutorial is a good starting point if you want to understand how HistFactory +## works. +## +## \date November 2026 +## \author Jonas Rembser, CERN + +import math + +import ROOT + + +def create_constraints(gammas, samples, global_observables): + constraints = ROOT.RooArgSet() + for i, gamma in enumerate(gammas): + nom_name = "nom_" + gamma.GetName() + tau_name = gamma.GetName() + "_tau" + mean_name = gamma.GetName() + "_mean" + poisson_name = gamma.GetName() + "_constraint" + + nom_val = sum(func.dataHist().weight(i) for func in samples) + + # Like in HistFactory, the gamma of an empty bin is fixed and not + # constrained. Otherwise, it would be a flat direction of the likelihood. + # The other gammas get the relative Poisson uncertainty as initial error. + if nom_val == 0.0: + gamma.setConstant(True) + continue + gamma.setError(1.0 / math.sqrt(nom_val)) + + nominal = ROOT.RooRealVar(nom_name, nom_name, nom_val, 0.0, ROOT.RooNumber.infinity()) + nominal.setConstant(True) + global_observables.add(nominal) + + tau = ROOT.RooConstVar(tau_name, tau_name, nominal.getVal()) + + mean = ROOT.RooProduct(mean_name, mean_name, gamma, tau) + + poisson = ROOT.RooPoisson(poisson_name, poisson_name, nominal, mean) + + # The constraint pdf takes ownership of the objects it depends on, and + # the returned set takes ownership of the constraint pdf. + for helper in (nominal, tau, mean): + ROOT.SetOwnership(helper, False) + poisson.addOwnedComponents(ROOT.RooArgSet(nominal, tau, mean)) + + constraints.addOwned(poisson) + return constraints + + +ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.NumIntegration) + +# First, construct a likelihood model with a Gaussian signal on top of a uniform background +x = ROOT.RooRealVar("x", "x", -20, 20) +x.setBins(25) + +meanG = ROOT.RooRealVar("meanG", "meanG", 1, -10, 10) +sigG = ROOT.RooRealVar("sigG", "sigG", 1.5, -10, 10) +g = ROOT.RooGaussian("g", "Gauss", x, meanG, sigG) +u = ROOT.RooUniform("u", "Uniform", x) + +# Generate the data to be fitted +sigData = g.generate(x, 50) +bkgData = u.generate(x, 1000) + +sumData = ROOT.RooDataSet("sumData", "Gauss + Uniform", x) +sumData.append(sigData) +sumData.append(bkgData) + +# Make histogram templates for signal and background. +# Let's take a signal distribution with low statistics and a more accurate +# background distribution. +# Normally, these come from Monte Carlo simulations, but we will just generate them. +dh_sig = g.generateBinned(x, 50) +dh_bkg = u.generateBinned(x, 10000) + +# ***** Case 0 - 'Rigid templates' ***** + +# Construct histogram shapes for signal and background +p_h_sig = ROOT.RooHistFunc("p_h_sig", "p_h_sig", x, dh_sig) +p_h_bkg = ROOT.RooHistFunc("p_h_bkg", "p_h_bkg", x, dh_bkg) + +# Construct scale factors for adding the two distributions +Asig0 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) +Abkg0 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) + +# Construct the sum model +model0 = ROOT.RooRealSumPdf("model0", "model0", [p_h_sig, p_h_bkg], [Asig0, Abkg0], True) + +# ***** Case 1 - 'Barlow Beeston' ***** + +# Construct parameterized histogram shapes for signal and background +gammas_sig1 = ROOT.ParamHistFunc.createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10) +gammas_bkg1 = ROOT.ParamHistFunc.createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10) + +p_phf_sig1 = ROOT.ParamHistFunc("p_phf_sig1", "p_phf_sig1", [x], gammas_sig1) +p_phf_bkg1 = ROOT.ParamHistFunc("p_phf_bkg1", "p_phf_bkg1", [x], gammas_bkg1) + +p_p_sig1 = ROOT.RooProduct("p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1) +p_p_bkg1 = ROOT.RooProduct("p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1) + +Asig1 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) +Abkg1 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) + +# Construct the sum of these +model1_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig1, p_p_bkg1], [Asig1, Abkg1], True) + +# Construct the subsidiary poisson measurements constraining the histogram parameters +# These ensure that the bin contents of the histograms are only allowed to vary within +# the statistical uncertainty of the Monte Carlo. +globs1 = ROOT.RooArgSet() +constraints_sig1 = create_constraints(gammas_sig1, [p_h_sig], globs1) +constraints_bkg1 = create_constraints(gammas_bkg1, [p_h_bkg], globs1) + +# Construct the joint model with template PDFs and constraints +model1 = ROOT.RooProdPdf("model1", "model1", [model1_tmp, *constraints_sig1, *constraints_bkg1]) + +# ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** + +# Construct the histogram shapes, using the same parameters for signal and background. +# There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to +# fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background +# histogram, so the bins fluctuate synchronously. This reduces the number of parameters. +gammas2 = ROOT.ParamHistFunc.createParamSet("p_ph_gamma", x.numBins(), 0, 10) + +p_phf2 = ROOT.ParamHistFunc("p_phf2", "p_phf2", [x], gammas2) + +p_p_sig2 = ROOT.RooProduct("p_p_sig2", "p_p_sig2", p_h_sig, p_phf2) +p_p_bkg2 = ROOT.RooProduct("p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2) + +Asig2 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) +Abkg2 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) + +# As before, construct the sum of signal2 and background2 +model2_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig2, p_p_bkg2], [Asig2, Abkg2], True) + +# Construct the subsidiary poisson measurements constraining the statistical fluctuations +globs2 = ROOT.RooArgSet() +constraints2 = create_constraints(gammas2, [p_h_sig, p_h_bkg], globs2) + +# Construct the joint model +model2 = ROOT.RooProdPdf("model2", "model2", [model2_tmp, *constraints2]) + +# ************ Fit all models to data and plot ********************* + +result0 = model0.fitTo(sumData, PrintLevel=0, Save=True) +result1 = model1.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs1) +result2 = model2.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs2) + +can = ROOT.TCanvas("can", "", 1500, 600) +can.Divide(3, 1) + +pt = ROOT.TPaveText(-19.5, 1, -2, 25) +pt.SetFillStyle(0) +pt.SetBorderSize(0) + +can.cd(1) +frame = x.frame(Title="No template uncertainties") +# Plot data to enable automatic determination of model0 normalisation: +sumData.plotOn(frame) +model0.plotOn(frame, LineColor="b", VisualizeError=result0) +# Plot data again to show it on top of model0 error bands: +sumData.plotOn(frame) +# Plot model components +model0.plotOn(frame, LineColor="b") +model0.plotOn(frame, Components=[p_h_sig], LineColor="kAzure") +model0.plotOn(frame, Components=[p_h_bkg], LineColor="r") +model0.paramOn(frame) + +sigData.plotOn(frame, MarkerColor="b") +frame.Draw() + +for text in [ + "No template uncertainties", + "are taken into account.", + "This leads to low errors", + "for the parameters A, since", + "the only source of errors", + "are the data statistics.", +]: + pt.AddText(text) +pt.DrawClone() + +can.cd(2) +frame = x.frame(Title="Barlow Beeston for Sig & Bkg separately") +sumData.plotOn(frame) +model1.plotOn(frame, LineColor="b", VisualizeError=result1) +# Plot data again to show it on top of error bands: +sumData.plotOn(frame) +model1.plotOn(frame, LineColor="b") +model1.plotOn(frame, Components=[p_p_sig1], LineColor="kAzure") +model1.plotOn(frame, Components=[p_p_bkg1], LineColor="r") +model1.paramOn(frame, Parameters=[Asig1, Abkg1]) + +sigData.plotOn(frame, MarkerColor="b") +frame.Draw() + +pt.Clear() +for text in [ + "With gamma parameters, the", + "signal & background templates", + "can adapt to the data.", + "Note how the blue signal", + "template changes its shape.", + "This leads to higher errors", + "of the scale parameters A.", +]: + pt.AddText(text) +pt.DrawClone() + +can.cd(3) +frame = x.frame(Title="Barlow Beeston light for (Sig+Bkg)") +sumData.plotOn(frame) +model2.plotOn(frame, LineColor="b", VisualizeError=result2) +# Plot data again to show it on top of model0 error bands: +sumData.plotOn(frame) +model2.plotOn(frame, LineColor="b") +model2.plotOn(frame, Components=[p_p_sig2], LineColor="kAzure") +model2.plotOn(frame, Components=[p_p_bkg2], LineColor="r") +model2.paramOn(frame, Parameters=[Asig2, Abkg2]) + +sigData.plotOn(frame, MarkerColor="b") +frame.Draw() + +pt.Clear() +for text in [ + "When signal and background", + "template share one gamma para-", + "meter per bin, they adapt less.", + "The errors of the A parameters", + "also shrink slightly.", +]: + pt.AddText(text) +pt.DrawClone() + +print(f"Asig [normal ] = {Asig0.getVal()} +/- {Asig0.getError()}") +print(f"Asig [BB ] = {Asig1.getVal()} +/- {Asig1.getError()}") +print(f"Asig [BBlight] = {Asig2.getVal()} +/- {Asig2.getError()}") + +can.SaveAs("rf709b_BarlowBeestonHistFactory.png") From f478d9d1edd796d399af2f501c8fd565ecfd0189 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 6 Nov 2026 11:08:23 +0000 Subject: [PATCH 4/5] [RF] Replace rf709 tutorials with the HistFactory-style rewrites Model the template uncertainties with ParamHistFunc and per-bin RooPoisson constraints, like HistFactory, instead of RooParamHistFunc and RooHistConstraint, and remove the rf709b copies. The docs of these two classes no longer point to rf709, since it doesn't use them anymore. Assisted-by: claude-code:claude-opus-5-5 --- roofit/roofit/src/RooHistConstraint.cxx | 2 - roofit/roofit/src/RooParamHistFunc.cxx | 2 - tutorials/roofit/roofit/index.md | 1 - tutorials/roofit/roofit/rf709_BarlowBeeston.C | 404 +++++++++--------- .../roofit/roofit/rf709_BarlowBeeston.py | 139 +++--- .../roofit/rf709b_BarlowBeestonHistFactory.C | 263 ------------ .../roofit/rf709b_BarlowBeestonHistFactory.py | 251 ----------- 7 files changed, 300 insertions(+), 762 deletions(-) delete mode 100644 tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.C delete mode 100644 tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py diff --git a/roofit/roofit/src/RooHistConstraint.cxx b/roofit/roofit/src/RooHistConstraint.cxx index 80cfc3f7113d9..039037383f141 100644 --- a/roofit/roofit/src/RooHistConstraint.cxx +++ b/roofit/roofit/src/RooHistConstraint.cxx @@ -16,8 +16,6 @@ * * It can therefore be used to estimate the Monte Carlo uncertainty of a fit. * - * Check also the tutorial rf709_BarlowBeeston.C - * */ #include diff --git a/roofit/roofit/src/RooParamHistFunc.cxx b/roofit/roofit/src/RooParamHistFunc.cxx index 3f37cd548ece1..778ce475bdee3 100644 --- a/roofit/roofit/src/RooParamHistFunc.cxx +++ b/roofit/roofit/src/RooParamHistFunc.cxx @@ -19,8 +19,6 @@ * The \f$ \gamma_i \f$ can therefore be used to parametrise statistical uncertainties of the histogram * template. In conjunction with a constraint term, this can be used to implement the Barlow-Beeston method. * The constraint can be implemented using RooHistConstraint. - * - * See also the tutorial rf709_BarlowBeeston.C */ #include diff --git a/tutorials/roofit/roofit/index.md b/tutorials/roofit/roofit/index.md index f4ea1335e6558..347248a7f765f 100644 --- a/tutorials/roofit/roofit/index.md +++ b/tutorials/roofit/roofit/index.md @@ -180,7 +180,6 @@ Explore the tutorials below to discover the main features of RooFit. A more inde | | rf617_simulation_based_inference_multidimensional.py | Use Simulation Based Inference (SBI) in multiple dimensions in RooFit.| | | rf618_mixture_models.py | Use of mixture models in RooFit.| | rf709_BarlowBeeston.C| rf709_BarlowBeeston.py | Implementing the Barlow-Beeston method for taking into account the statistical uncertainty of a Monte-Carlo fit template.| -| rf709b_BarlowBeestonHistFactory.C| rf709b_BarlowBeestonHistFactory.py | Implementing the Barlow-Beeston method with the same kind of RooFit models that HistFactory builds.| | rf710_roopoly.C| rf710_roopoly.py | Taylor expansion of RooFit functions using the taylorExpand function with RooPolyFunc.| | rf711_lagrangianmorph.C| rf711_lagrangianmorph.py | Morphing effective field theory distributions with RooLagrangianMorphFunc. A morphing function as a function of one coefficient is setup and can be used to obtain the distribution for any value of the coefficient.| | rf712_lagrangianmorphfit.C| rf712_lagrangianmorphfit.py | Performing a simple fit with RooLagrangianMorphFunc. A morphing function is setup as a function of three variables and a fit is performed on a pseudo-dataset.| diff --git a/tutorials/roofit/roofit/rf709_BarlowBeeston.C b/tutorials/roofit/roofit/rf709_BarlowBeeston.C index f472e3f0da92d..100badec70329 100644 --- a/tutorials/roofit/roofit/rf709_BarlowBeeston.C +++ b/tutorials/roofit/roofit/rf709_BarlowBeeston.C @@ -11,7 +11,7 @@ /// Based on a demo by Wouter Verkerke // /// \date June 2019 -/// \author Stephan Hageboeck, CERN +/// \author Jonas Rembser, Stephan Hageboeck, CERN #include "RooRealVar.h" #include "RooGaussian.h" @@ -20,219 +20,241 @@ #include "RooDataHist.h" #include "RooHistFunc.h" #include "RooRealSumPdf.h" -#include "RooParamHistFunc.h" -#include "RooHistConstraint.h" +#include "RooProduct.h" +#include "RooConstVar.h" +#include "RooPoisson.h" #include "RooProdPdf.h" +#include "RooFitResult.h" #include "RooPlot.h" +#include "RooMsgService.h" +#include "RooStats/HistFactory/ParamHistFunc.h" #include "TCanvas.h" #include "TPaveText.h" +#include #include #include -using namespace RooFit; +RooArgSet createConstraints(RooArgList const &gammas, RooArgList const &samples, RooArgSet &globalObservables) +{ + RooArgSet constraints; + for (std::size_t i = 0; i < gammas.size(); ++i) { + auto *gamma = static_cast(&gammas[i]); + + auto nomName = "nom_" + std::string(gamma->GetName()); + auto tauName = std::string(gamma->GetName()) + "_tau"; + auto meanName = std::string(gamma->GetName()) + "_mean"; + auto poissonName = std::string(gamma->GetName()) + "_constraint"; + + double nomVal = 0.0; + for (auto *func : static_range_cast(samples)) { + nomVal += func->dataHist().weight(i); + } + + // Like in HistFactory, the gamma of an empty bin is fixed and not + // constrained. Otherwise, it would be a flat direction of the likelihood. + // The other gammas get the relative Poisson uncertainty as initial error. + if (nomVal == 0.0) { + gamma->setConstant(true); + continue; + } + gamma->setError(1.0 / std::sqrt(nomVal)); + + auto nominal = std::make_unique(nomName.c_str(), nomName.c_str(), nomVal, 0.0, RooNumber::infinity()); + nominal->setConstant(true); + globalObservables.add(*nominal); + + auto tau = std::make_unique(tauName.c_str(), tauName.c_str(), nominal->getVal()); + + auto mean = std::make_unique(meanName.c_str(), meanName.c_str(), *gamma, *tau); + mean->addOwnedComponents(std::move(tau)); + + auto poisson = std::make_unique(poissonName.c_str(), poissonName.c_str(), *nominal, *mean); + poisson->addOwnedComponents(std::move(nominal)); + poisson->addOwnedComponents(std::move(mean)); + + constraints.addOwned(std::move(poisson)); + } + return constraints; +} void rf709_BarlowBeeston() { - // First, construct a likelihood model with a Gaussian signal on top of a uniform background - RooRealVar x("x", "x", -20, 20); - x.setBins(25); - - RooRealVar meanG("meanG", "meanG", 1, -10, 10); - RooRealVar sigG("sigG", "sigG", 1.5, -10, 10); - RooGaussian g("g", "Gauss", x, meanG, sigG); - RooUniform u("u", "Uniform", x); - - - // Generate the data to be fitted - std::unique_ptr sigData(g.generate(x, 50)); - std::unique_ptr bkgData(u.generate(x, 1000)); - - RooDataSet sumData("sumData", "Gauss + Uniform", x); - sumData.append(*sigData); - sumData.append(*bkgData); - - - // Make histogram templates for signal and background. - // Let's take a signal distribution with low statistics and a more accurate - // background distribution. - // Normally, these come from Monte Carlo simulations, but we will just generate them. - std::unique_ptr dh_sig( g.generateBinned(x, 50) ); - std::unique_ptr dh_bkg( u.generateBinned(x, 10000) ); - - - // ***** Case 0 - 'Rigid templates' ***** - - // Construct histogram shapes for signal and background - RooHistFunc p_h_sig("p_h_sig","p_h_sig",x,*dh_sig); - RooHistFunc p_h_bkg("p_h_bkg","p_h_bkg",x,*dh_bkg); - - // Construct scale factors for adding the two distributions - RooRealVar Asig0("Asig","Asig",1,0.01,5000); - RooRealVar Abkg0("Abkg","Abkg",1,0.01,5000); - - // Construct the sum model - RooRealSumPdf model0("model0","model0", - RooArgList(p_h_sig,p_h_bkg), - RooArgList(Asig0,Abkg0), - true); + using namespace RooFit; + RooMsgService::instance().getStream(1).removeTopic(NumIntegration); + // First, construct a likelihood model with a Gaussian signal on top of a uniform background + RooRealVar x("x", "x", -20, 20); + x.setBins(25); - // ***** Case 1 - 'Barlow Beeston' ***** + RooRealVar meanG("meanG", "meanG", 1, -10, 10); + RooRealVar sigG("sigG", "sigG", 1.5, -10, 10); + RooGaussian g("g", "Gauss", x, meanG, sigG); + RooUniform u("u", "Uniform", x); - // Construct parameterized histogram shapes for signal and background - RooParamHistFunc p_ph_sig1("p_ph_sig","p_ph_sig",*dh_sig, x); - RooParamHistFunc p_ph_bkg1("p_ph_bkg","p_ph_bkg",*dh_bkg, x); + // Generate the data to be fitted + std::unique_ptr sigData(g.generate(x, 50)); + std::unique_ptr bkgData(u.generate(x, 1000)); - RooRealVar Asig1("Asig","Asig",1,0.01,5000); - RooRealVar Abkg1("Abkg","Abkg",1,0.01,5000); + RooDataSet sumData("sumData", "Gauss + Uniform", x); + sumData.append(*sigData); + sumData.append(*bkgData); - // Construct the sum of these - RooRealSumPdf model_tmp("sp_ph", "sp_ph", - RooArgList(p_ph_sig1,p_ph_bkg1), - RooArgList(Asig1,Abkg1), - true); + // Make histogram templates for signal and background. + // Let's take a signal distribution with low statistics and a more accurate + // background distribution. + // Normally, these come from Monte Carlo simulations, but we will just generate them. + std::unique_ptr dh_sig(g.generateBinned(x, 50)); + std::unique_ptr dh_bkg(u.generateBinned(x, 10000)); - // Construct the subsidiary poisson measurements constraining the histogram parameters - // These ensure that the bin contents of the histograms are only allowed to vary within - // the statistical uncertainty of the Monte Carlo. - RooHistConstraint hc_sig("hc_sig","hc_sig",p_ph_sig1); - RooHistConstraint hc_bkg("hc_bkg","hc_bkg",p_ph_bkg1); + // ***** Case 0 - 'Rigid templates' ***** - // Construct the joint model with template PDFs and constraints - RooProdPdf model1("model1","model1",RooArgSet(hc_sig,hc_bkg),Conditional(model_tmp,x)); + // Construct histogram shapes for signal and background + RooHistFunc p_h_sig("p_h_sig", "p_h_sig", x, *dh_sig); + RooHistFunc p_h_bkg("p_h_bkg", "p_h_bkg", x, *dh_bkg); + // Construct scale factors for adding the two distributions + RooRealVar Asig0("Asig", "Asig", 1, 0.01, 5000); + RooRealVar Abkg0("Abkg", "Abkg", 1, 0.01, 5000); + // Construct the sum model + RooRealSumPdf model0("model0", "model0", {p_h_sig, p_h_bkg}, {Asig0, Abkg0}, true); - // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** + // ***** Case 1 - 'Barlow Beeston' ***** - // Construct the histogram shapes, using the same parameters for signal and background - // This requires passing the first histogram to the second, so that their common parameters - // can be re-used. - // The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`. - // This allows bin 0 to fluctuate up and down. - // Then, the SAME parameters are connected to the background histogram, so the bins fluctuate - // synchronously. This reduces the number of parameters. - RooParamHistFunc p_ph_sig2("p_ph_sig2", "p_ph_sig2", *dh_sig, x); - RooParamHistFunc p_ph_bkg2("p_ph_bkg2", "p_ph_bkg2", *dh_bkg, x, &p_ph_sig2, true); - - RooRealVar Asig2("Asig","Asig",1,0.01,5000); - RooRealVar Abkg2("Abkg","Abkg",1,0.01,5000); - - // As before, construct the sum of signal2 and background2 - RooRealSumPdf model2_tmp("sp_ph","sp_ph", - RooArgList(p_ph_sig2,p_ph_bkg2), - RooArgList(Asig2,Abkg2), - true); - - // Construct the subsidiary poisson measurements constraining the statistical fluctuations - RooHistConstraint hc_sigbkg("hc_sigbkg","hc_sigbkg",RooArgSet(p_ph_sig2,p_ph_bkg2)); - - // Construct the joint model - RooProdPdf model2("model2","model2",hc_sigbkg, RooFit::Conditional(model2_tmp,x)); - - - - // ************ Fit all models to data and plot ********************* - - auto result0 = model0.fitTo(sumData, PrintLevel(0), Save()); - auto result1 = model1.fitTo(sumData, PrintLevel(0), Save()); - auto result2 = model2.fitTo(sumData, PrintLevel(0), Save()); - - - TCanvas* can = new TCanvas("can", "", 1500, 600); - can->Divide(3,1); - - TPaveText pt(-19.5, 1, -2, 25); - pt.SetFillStyle(0); - pt.SetBorderSize(0); - - - can->cd(1); - auto frame = x.frame(Title("No template uncertainties")); - // Plot data to enable automatic determination of model0 normalisation: - sumData.plotOn(frame); - model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0)); - // Plot data again to show it on top of model0 error bands: - sumData.plotOn(frame); - // Plot model components - model0.plotOn(frame, LineColor(kBlue)); - model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure)); - model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed)); - model0.paramOn(frame); - - sigData->plotOn(frame, MarkerColor(kBlue)); - frame->Draw(); - - for (auto text : { - "No template uncertainties", - "are taken into account.", - "This leads to low errors", - "for the parameters A, since", - "the only source of errors", - "are the data statistics."}) { - pt.AddText(text); - } - pt.DrawClone(); - - - can->cd(2); - frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately")); - sumData.plotOn(frame); - model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1)); - // Plot data again to show it on top of error bands: - sumData.plotOn(frame); - model1.plotOn(frame, LineColor(kBlue)); - model1.plotOn(frame, Components(p_ph_sig1), LineColor(kAzure)); - model1.plotOn(frame, Components(p_ph_bkg1), LineColor(kRed)); - model1.paramOn(frame, Parameters(RooArgSet(Asig1, Abkg1))); - - sigData->plotOn(frame, MarkerColor(kBlue)); - frame->Draw(); - - pt.Clear(); - for (auto text : { - "With gamma parameters, the", - "signal & background templates", - "can adapt to the data.", - "Note how the blue signal", - "template changes its shape.", - "This leads to higher errors", - "of the scale parameters A."}) { - pt.AddText(text); - } - pt.DrawClone(); - - can->cd(3); - frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)")); - sumData.plotOn(frame); - model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2)); - // Plot data again to show it on top of model0 error bands: - sumData.plotOn(frame); - model2.plotOn(frame, LineColor(kBlue)); - model2.plotOn(frame, Components(p_ph_sig2), LineColor(kAzure)); - model2.plotOn(frame, Components(p_ph_bkg2), LineColor(kRed)); - model2.paramOn(frame, Parameters(RooArgSet(Asig2, Abkg2))); - - sigData->plotOn(frame, MarkerColor(kBlue)); - frame->Draw(); - - pt.Clear(); - for (auto text : { - "When signal and background", - "template share one gamma para-", - "meter per bin, they adapt less.", - "The errors of the A parameters", - "also shrink slightly."}) { - pt.AddText(text); - } - pt.DrawClone(); - - - std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl; - std::cout << "Asig [BB ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl; - std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl; + // Construct parameterized histogram shapes for signal and background + RooArgList gammas_sig1{ParamHistFunc::createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10)}; + RooArgList gammas_bkg1{ParamHistFunc::createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10)}; + ParamHistFunc p_phf_sig1{"p_phf_sig1", "p_phf_sig1", x, gammas_sig1}; + ParamHistFunc p_phf_bkg1{"p_phf_bkg1", "p_phf_bkg1", x, gammas_bkg1}; + + RooProduct p_p_sig1{"p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1}; + RooProduct p_p_bkg1{"p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1}; + + RooRealVar Asig1("Asig", "Asig", 1, 0.01, 5000); + RooRealVar Abkg1("Abkg", "Abkg", 1, 0.01, 5000); + + // Construct the sum of these + RooRealSumPdf model1_tmp("sp_ph", "sp_ph", {p_p_sig1, p_p_bkg1}, {Asig1, Abkg1}, true); + + // Construct the subsidiary poisson measurements constraining the histogram parameters + // These ensure that the bin contents of the histograms are only allowed to vary within + // the statistical uncertainty of the Monte Carlo. + RooArgSet globs1; + RooArgSet constraints_sig1{createConstraints(gammas_sig1, p_h_sig, globs1)}; + RooArgSet constraints_bkg1{createConstraints(gammas_bkg1, p_h_bkg, globs1)}; + + // Construct the joint model with template PDFs and constraints + RooProdPdf model1("model1", "model1", {model1_tmp, constraints_sig1, constraints_bkg1}); + + // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** + + // Construct the histogram shapes, using the same parameters for signal and background. + // There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to + // fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background + // histogram, so the bins fluctuate synchronously. This reduces the number of parameters. + RooArgList gammas2{ParamHistFunc::createParamSet("p_ph_gamma", x.numBins(), 0, 10)}; + + ParamHistFunc p_phf2{"p_phf2", "p_phf2", x, gammas2}; + + RooProduct p_p_sig2{"p_p_sig2", "p_p_sig2", p_h_sig, p_phf2}; + RooProduct p_p_bkg2{"p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2}; + + RooRealVar Asig2("Asig", "Asig", 1, 0.01, 5000); + RooRealVar Abkg2("Abkg", "Abkg", 1, 0.01, 5000); + + // As before, construct the sum of signal2 and background2 + RooRealSumPdf model2_tmp("sp_ph", "sp_ph", {p_p_sig2, p_p_bkg2}, {Asig2, Abkg2}, true); + + // Construct the subsidiary poisson measurements constraining the statistical fluctuations + RooArgSet globs2; + RooArgSet constraints2{createConstraints(gammas2, {p_h_sig, p_h_bkg}, globs2)}; + + // Construct the joint model + RooProdPdf model2("model2", "model2", {model2_tmp, constraints2}); + + // ************ Fit all models to data and plot ********************* + + std::unique_ptr result0{model0.fitTo(sumData, PrintLevel(0), Save())}; + std::unique_ptr result1{model1.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs1))}; + std::unique_ptr result2{model2.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs2))}; + + TCanvas *can = new TCanvas("can", "", 1500, 600); + can->Divide(3, 1); + + TPaveText pt(-19.5, 1, -2, 25); + pt.SetFillStyle(0); + pt.SetBorderSize(0); + + can->cd(1); + auto frame = x.frame(Title("No template uncertainties")); + // Plot data to enable automatic determination of model0 normalisation: + sumData.plotOn(frame); + model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0)); + // Plot data again to show it on top of model0 error bands: + sumData.plotOn(frame); + // Plot model components + model0.plotOn(frame, LineColor(kBlue)); + model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure)); + model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed)); + model0.paramOn(frame); + + sigData->plotOn(frame, MarkerColor(kBlue)); + frame->Draw(); + + for (auto text : {"No template uncertainties", "are taken into account.", "This leads to low errors", + "for the parameters A, since", "the only source of errors", "are the data statistics."}) { + pt.AddText(text); + } + pt.DrawClone(); + + can->cd(2); + frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately")); + sumData.plotOn(frame); + model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1)); + // Plot data again to show it on top of error bands: + sumData.plotOn(frame); + model1.plotOn(frame, LineColor(kBlue)); + model1.plotOn(frame, Components(p_p_sig1), LineColor(kAzure)); + model1.plotOn(frame, Components(p_p_bkg1), LineColor(kRed)); + model1.paramOn(frame, Parameters({Asig1, Abkg1})); + + sigData->plotOn(frame, MarkerColor(kBlue)); + frame->Draw(); + + pt.Clear(); + for (auto text : {"With gamma parameters, the", "signal & background templates", "can adapt to the data.", + "Note how the blue signal", "template changes its shape.", "This leads to higher errors", + "of the scale parameters A."}) { + pt.AddText(text); + } + pt.DrawClone(); + + can->cd(3); + frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)")); + sumData.plotOn(frame); + model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2)); + // Plot data again to show it on top of model0 error bands: + sumData.plotOn(frame); + model2.plotOn(frame, LineColor(kBlue)); + model2.plotOn(frame, Components(p_p_sig2), LineColor(kAzure)); + model2.plotOn(frame, Components(p_p_bkg2), LineColor(kRed)); + model2.paramOn(frame, Parameters({Asig2, Abkg2})); + + sigData->plotOn(frame, MarkerColor(kBlue)); + frame->Draw(); + + pt.Clear(); + for (auto text : {"When signal and background", "template share one gamma para-", "meter per bin, they adapt less.", + "The errors of the A parameters", "also shrink slightly."}) { + pt.AddText(text); + } + pt.DrawClone(); + + std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl; + std::cout << "Asig [BB ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl; + std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl; } diff --git a/tutorials/roofit/roofit/rf709_BarlowBeeston.py b/tutorials/roofit/roofit/rf709_BarlowBeeston.py index fff319beadefe..0151c8e3e2b63 100644 --- a/tutorials/roofit/roofit/rf709_BarlowBeeston.py +++ b/tutorials/roofit/roofit/rf709_BarlowBeeston.py @@ -10,11 +10,53 @@ ## ## Based on a demo by Wouter Verkerke ## \date June 2021 -## \author Harshal Shende, Stephan Hageboeck (C++ version) +## \author Jonas Rembser, Harshal Shende, Stephan Hageboeck (C++ version) +import math import ROOT + +def create_constraints(gammas, samples, global_observables): + constraints = ROOT.RooArgSet() + for i, gamma in enumerate(gammas): + nom_name = "nom_" + gamma.GetName() + tau_name = gamma.GetName() + "_tau" + mean_name = gamma.GetName() + "_mean" + poisson_name = gamma.GetName() + "_constraint" + + nom_val = sum(func.dataHist().weight(i) for func in samples) + + # Like in HistFactory, the gamma of an empty bin is fixed and not + # constrained. Otherwise, it would be a flat direction of the likelihood. + # The other gammas get the relative Poisson uncertainty as initial error. + if nom_val == 0.0: + gamma.setConstant(True) + continue + gamma.setError(1.0 / math.sqrt(nom_val)) + + nominal = ROOT.RooRealVar(nom_name, nom_name, nom_val, 0.0, ROOT.RooNumber.infinity()) + nominal.setConstant(True) + global_observables.add(nominal) + + tau = ROOT.RooConstVar(tau_name, tau_name, nominal.getVal()) + + mean = ROOT.RooProduct(mean_name, mean_name, gamma, tau) + + poisson = ROOT.RooPoisson(poisson_name, poisson_name, nominal, mean) + + # The constraint pdf takes ownership of the objects it depends on, and + # the returned set takes ownership of the constraint pdf. + for helper in (nominal, tau, mean): + ROOT.SetOwnership(helper, False) + poisson.addOwnedComponents(ROOT.RooArgSet(nominal, tau, mean)) + + constraints.addOwned(poisson) + return constraints + + +ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.NumIntegration) + # First, construct a likelihood model with a Gaussian signal on top of a uniform background x = ROOT.RooRealVar("x", "x", -20, 20) x.setBins(25) @@ -24,7 +66,6 @@ g = ROOT.RooGaussian("g", "Gauss", x, meanG, sigG) u = ROOT.RooUniform("u", "Uniform", x) - # Generate the data to be fitted sigData = g.generate(x, 50) bkgData = u.generate(x, 1000) @@ -33,7 +74,6 @@ sumData.append(sigData) sumData.append(bkgData) - # Make histogram templates for signal and background. # Let's take a signal distribution with low statistics and a more accurate # background distribution. @@ -41,8 +81,7 @@ dh_sig = g.generateBinned(x, 50) dh_bkg = u.generateBinned(x, 10000) - -# Case 0 - 'Rigid templates' +# ***** Case 0 - 'Rigid templates' ***** # Construct histogram shapes for signal and background p_h_sig = ROOT.RooHistFunc("p_h_sig", "p_h_sig", x, dh_sig) @@ -55,60 +94,65 @@ # Construct the sum model model0 = ROOT.RooRealSumPdf("model0", "model0", [p_h_sig, p_h_bkg], [Asig0, Abkg0], True) - -# Case 1 - 'Barlow Beeston' +# ***** Case 1 - 'Barlow Beeston' ***** # Construct parameterized histogram shapes for signal and background -p_ph_sig1 = ROOT.RooParamHistFunc("p_ph_sig", "p_ph_sig", dh_sig, x) -p_ph_bkg1 = ROOT.RooParamHistFunc("p_ph_bkg", "p_ph_bkg", dh_bkg, x) +gammas_sig1 = ROOT.ParamHistFunc.createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10) +gammas_bkg1 = ROOT.ParamHistFunc.createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10) + +p_phf_sig1 = ROOT.ParamHistFunc("p_phf_sig1", "p_phf_sig1", [x], gammas_sig1) +p_phf_bkg1 = ROOT.ParamHistFunc("p_phf_bkg1", "p_phf_bkg1", [x], gammas_bkg1) + +p_p_sig1 = ROOT.RooProduct("p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1) +p_p_bkg1 = ROOT.RooProduct("p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1) Asig1 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) Abkg1 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) # Construct the sum of these -model_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_ph_sig1, p_ph_bkg1], [Asig1, Abkg1], True) +model1_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig1, p_p_bkg1], [Asig1, Abkg1], True) # Construct the subsidiary poisson measurements constraining the histogram parameters # These ensure that the bin contents of the histograms are only allowed to vary within # the statistical uncertainty of the Monte Carlo. -hc_sig = ROOT.RooHistConstraint("hc_sig", "hc_sig", p_ph_sig1) -hc_bkg = ROOT.RooHistConstraint("hc_bkg", "hc_bkg", p_ph_bkg1) +globs1 = ROOT.RooArgSet() +constraints_sig1 = create_constraints(gammas_sig1, [p_h_sig], globs1) +constraints_bkg1 = create_constraints(gammas_bkg1, [p_h_bkg], globs1) # Construct the joint model with template PDFs and constraints -model1 = ROOT.RooProdPdf("model1", "model1", [hc_sig, hc_bkg], Conditional=(model_tmp, x)) +model1 = ROOT.RooProdPdf("model1", "model1", [model1_tmp, *constraints_sig1, *constraints_bkg1]) +# ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** -# Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) +# Construct the histogram shapes, using the same parameters for signal and background. +# There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to +# fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background +# histogram, so the bins fluctuate synchronously. This reduces the number of parameters. +gammas2 = ROOT.ParamHistFunc.createParamSet("p_ph_gamma", x.numBins(), 0, 10) -# Construct the histogram shapes, using the same parameters for signal and background -# This requires passing the first histogram to the second, so that their common parameters -# can be re-used. -# The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`. -# This allows bin 0 to fluctuate up and down. -# Then, the SAME parameters are connected to the background histogram, so the bins fluctuate -# synchronously. This reduces the number of parameters. -p_ph_sig2 = ROOT.RooParamHistFunc("p_ph_sig2", "p_ph_sig2", dh_sig, x) -p_ph_bkg2 = ROOT.RooParamHistFunc("p_ph_bkg2", "p_ph_bkg2", dh_bkg, x, p_ph_sig2, True) +p_phf2 = ROOT.ParamHistFunc("p_phf2", "p_phf2", [x], gammas2) + +p_p_sig2 = ROOT.RooProduct("p_p_sig2", "p_p_sig2", p_h_sig, p_phf2) +p_p_bkg2 = ROOT.RooProduct("p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2) Asig2 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) Abkg2 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) # As before, construct the sum of signal2 and background2 -model2_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_ph_sig2, p_ph_bkg2], [Asig2, Abkg2], True) +model2_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig2, p_p_bkg2], [Asig2, Abkg2], True) # Construct the subsidiary poisson measurements constraining the statistical fluctuations -hc_sigbkg = ROOT.RooHistConstraint("hc_sigbkg", "hc_sigbkg", [p_ph_sig2, p_ph_bkg2]) +globs2 = ROOT.RooArgSet() +constraints2 = create_constraints(gammas2, [p_h_sig, p_h_bkg], globs2) # Construct the joint model -model2 = ROOT.RooProdPdf("model2", "model2", hc_sigbkg, Conditional=(model2_tmp, x)) - +model2 = ROOT.RooProdPdf("model2", "model2", [model2_tmp, *constraints2]) # ************ Fit all models to data and plot ********************* result0 = model0.fitTo(sumData, PrintLevel=0, Save=True) -result1 = model1.fitTo(sumData, PrintLevel=0, Save=True) -result2 = model2.fitTo(sumData, PrintLevel=0, Save=True) - +result1 = model1.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs1) +result2 = model2.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs2) can = ROOT.TCanvas("can", "", 1500, 600) can.Divide(3, 1) @@ -117,7 +161,6 @@ pt.SetFillStyle(0) pt.SetBorderSize(0) - can.cd(1) frame = x.frame(Title="No template uncertainties") # Plot data to enable automatic determination of model0 normalisation: @@ -134,20 +177,17 @@ sigData.plotOn(frame, MarkerColor="b") frame.Draw() -pt_text1 = [ +for text in [ "No template uncertainties", "are taken into account.", "This leads to low errors", "for the parameters A, since", "the only source of errors", "are the data statistics.", -] -for text in pt_text1: +]: pt.AddText(text) - pt.DrawClone() - can.cd(2) frame = x.frame(Title="Barlow Beeston for Sig & Bkg separately") sumData.plotOn(frame) @@ -155,15 +195,15 @@ # Plot data again to show it on top of error bands: sumData.plotOn(frame) model1.plotOn(frame, LineColor="b") -model1.plotOn(frame, Components=[p_ph_sig1], LineColor="kAzure") -model1.plotOn(frame, Components=[p_ph_bkg1], LineColor="r") +model1.plotOn(frame, Components=[p_p_sig1], LineColor="kAzure") +model1.plotOn(frame, Components=[p_p_bkg1], LineColor="r") model1.paramOn(frame, Parameters=[Asig1, Abkg1]) sigData.plotOn(frame, MarkerColor="b") frame.Draw() pt.Clear() -pt_text2 = [ +for text in [ "With gamma parameters, the", "signal & background templates", "can adapt to the data.", @@ -171,11 +211,8 @@ "template changes its shape.", "This leads to higher errors", "of the scale parameters A.", -] - -for text in pt_text2: +]: pt.AddText(text) - pt.DrawClone() can.cd(3) @@ -185,28 +222,26 @@ # Plot data again to show it on top of model0 error bands: sumData.plotOn(frame) model2.plotOn(frame, LineColor="b") -model2.plotOn(frame, Components=[p_ph_sig2], LineColor="kAzure") -model2.plotOn(frame, Components=[p_ph_bkg2], LineColor="r") +model2.plotOn(frame, Components=[p_p_sig2], LineColor="kAzure") +model2.plotOn(frame, Components=[p_p_bkg2], LineColor="r") model2.paramOn(frame, Parameters=[Asig2, Abkg2]) sigData.plotOn(frame, MarkerColor="b") frame.Draw() pt.Clear() -pt_text3 = [ +for text in [ "When signal and background", "template share one gamma para-", "meter per bin, they adapt less.", "The errors of the A parameters", "also shrink slightly.", -] -for text in pt_text3: +]: pt.AddText(text) pt.DrawClone() - -print("Asig [normal ] = {} +/- {}".format(Asig0.getVal(), Asig0.getError())) -print("Asig [BB ] = {} +/- {}".format(Asig1.getVal(), Asig1.getError())) -print("Asig [BBlight] = {} +/- {}".format(Asig2.getVal(), Asig2.getError())) +print(f"Asig [normal ] = {Asig0.getVal()} +/- {Asig0.getError()}") +print(f"Asig [BB ] = {Asig1.getVal()} +/- {Asig1.getError()}") +print(f"Asig [BBlight] = {Asig2.getVal()} +/- {Asig2.getError()}") can.SaveAs("rf709_BarlowBeeston.png") diff --git a/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.C b/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.C deleted file mode 100644 index 0c58cf5bb1eab..0000000000000 --- a/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.C +++ /dev/null @@ -1,263 +0,0 @@ -/// \file -/// \ingroup tutorial_roofit_main -/// \notebook -js -/// Implementing the Barlow-Beeston method for taking into account the statistical -/// uncertainty of a Monte-Carlo fit template. -/// -/// \macro_image -/// \macro_code -/// \macro_output -/// -/// New implementation of the original rf709_BarlowBeeston tutorial, using the -/// same kind of RooFit models that are built by HistFactory. Therefore, this -/// tutorial is a good starting point if you want to understand how HistFactory -/// works. -// -/// \date June 2023 -/// \author Jonas Rembser, CERN - -#include "RooRealVar.h" -#include "RooGaussian.h" -#include "RooUniform.h" -#include "RooDataSet.h" -#include "RooDataHist.h" -#include "RooHistFunc.h" -#include "RooRealSumPdf.h" -#include "RooProduct.h" -#include "RooConstVar.h" -#include "RooPoisson.h" -#include "RooProdPdf.h" -#include "RooFitResult.h" -#include "RooPlot.h" -#include "RooMsgService.h" -#include "RooStats/HistFactory/ParamHistFunc.h" - -#include "TCanvas.h" -#include "TPaveText.h" - -#include -#include -#include - -RooArgSet createConstraints(RooArgList const &gammas, RooArgList const &samples, RooArgSet &globalObservables) -{ - RooArgSet constraints; - for (std::size_t i = 0; i < gammas.size(); ++i) { - auto *gamma = static_cast(&gammas[i]); - - auto nomName = "nom_" + std::string(gamma->GetName()); - auto tauName = std::string(gamma->GetName()) + "_tau"; - auto meanName = std::string(gamma->GetName()) + "_mean"; - auto poissonName = std::string(gamma->GetName()) + "_constraint"; - - double nomVal = 0.0; - for (auto *func : static_range_cast(samples)) { - nomVal += func->dataHist().weight(i); - } - - // Like in HistFactory, the gamma of an empty bin is fixed and not - // constrained. Otherwise, it would be a flat direction of the likelihood. - // The other gammas get the relative Poisson uncertainty as initial error. - if (nomVal == 0.0) { - gamma->setConstant(true); - continue; - } - gamma->setError(1.0 / std::sqrt(nomVal)); - - auto nominal = std::make_unique(nomName.c_str(), nomName.c_str(), nomVal, 0.0, RooNumber::infinity()); - nominal->setConstant(true); - globalObservables.add(*nominal); - - auto tau = std::make_unique(tauName.c_str(), tauName.c_str(), nominal->getVal()); - - auto mean = std::make_unique(meanName.c_str(), meanName.c_str(), *gamma, *tau); - mean->addOwnedComponents(std::move(tau)); - - auto poisson = std::make_unique(poissonName.c_str(), poissonName.c_str(), *nominal, *mean); - poisson->addOwnedComponents(std::move(nominal)); - poisson->addOwnedComponents(std::move(mean)); - - constraints.addOwned(std::move(poisson)); - } - return constraints; -} - -void rf709b_BarlowBeestonHistFactory() -{ - using namespace RooFit; - - RooMsgService::instance().getStream(1).removeTopic(NumIntegration); - - // First, construct a likelihood model with a Gaussian signal on top of a uniform background - RooRealVar x("x", "x", -20, 20); - x.setBins(25); - - RooRealVar meanG("meanG", "meanG", 1, -10, 10); - RooRealVar sigG("sigG", "sigG", 1.5, -10, 10); - RooGaussian g("g", "Gauss", x, meanG, sigG); - RooUniform u("u", "Uniform", x); - - // Generate the data to be fitted - std::unique_ptr sigData(g.generate(x, 50)); - std::unique_ptr bkgData(u.generate(x, 1000)); - - RooDataSet sumData("sumData", "Gauss + Uniform", x); - sumData.append(*sigData); - sumData.append(*bkgData); - - // Make histogram templates for signal and background. - // Let's take a signal distribution with low statistics and a more accurate - // background distribution. - // Normally, these come from Monte Carlo simulations, but we will just generate them. - std::unique_ptr dh_sig(g.generateBinned(x, 50)); - std::unique_ptr dh_bkg(u.generateBinned(x, 10000)); - - // ***** Case 0 - 'Rigid templates' ***** - - // Construct histogram shapes for signal and background - RooHistFunc p_h_sig("p_h_sig", "p_h_sig", x, *dh_sig); - RooHistFunc p_h_bkg("p_h_bkg", "p_h_bkg", x, *dh_bkg); - - // Construct scale factors for adding the two distributions - RooRealVar Asig0("Asig", "Asig", 1, 0.01, 5000); - RooRealVar Abkg0("Abkg", "Abkg", 1, 0.01, 5000); - - // Construct the sum model - RooRealSumPdf model0("model0", "model0", {p_h_sig, p_h_bkg}, {Asig0, Abkg0}, true); - - // ***** Case 1 - 'Barlow Beeston' ***** - - // Construct parameterized histogram shapes for signal and background - RooArgList gammas_sig1{ParamHistFunc::createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10)}; - RooArgList gammas_bkg1{ParamHistFunc::createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10)}; - - ParamHistFunc p_phf_sig1{"p_phf_sig1", "p_phf_sig1", x, gammas_sig1}; - ParamHistFunc p_phf_bkg1{"p_phf_bkg1", "p_phf_bkg1", x, gammas_bkg1}; - - RooProduct p_p_sig1{"p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1}; - RooProduct p_p_bkg1{"p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1}; - - RooRealVar Asig1("Asig", "Asig", 1, 0.01, 5000); - RooRealVar Abkg1("Abkg", "Abkg", 1, 0.01, 5000); - - // Construct the sum of these - RooRealSumPdf model1_tmp("sp_ph", "sp_ph", {p_p_sig1, p_p_bkg1}, {Asig1, Abkg1}, true); - - // Construct the subsidiary poisson measurements constraining the histogram parameters - // These ensure that the bin contents of the histograms are only allowed to vary within - // the statistical uncertainty of the Monte Carlo. - RooArgSet globs1; - RooArgSet constraints_sig1{createConstraints(gammas_sig1, p_h_sig, globs1)}; - RooArgSet constraints_bkg1{createConstraints(gammas_bkg1, p_h_bkg, globs1)}; - - // Construct the joint model with template PDFs and constraints - RooProdPdf model1("model1", "model1", {model1_tmp, constraints_sig1, constraints_bkg1}); - - // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** - - // Construct the histogram shapes, using the same parameters for signal and background. - // There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to - // fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background - // histogram, so the bins fluctuate synchronously. This reduces the number of parameters. - RooArgList gammas2{ParamHistFunc::createParamSet("p_ph_gamma", x.numBins(), 0, 10)}; - - ParamHistFunc p_phf2{"p_phf2", "p_phf2", x, gammas2}; - - RooProduct p_p_sig2{"p_p_sig2", "p_p_sig2", p_h_sig, p_phf2}; - RooProduct p_p_bkg2{"p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2}; - - RooRealVar Asig2("Asig", "Asig", 1, 0.01, 5000); - RooRealVar Abkg2("Abkg", "Abkg", 1, 0.01, 5000); - - // As before, construct the sum of signal2 and background2 - RooRealSumPdf model2_tmp("sp_ph", "sp_ph", {p_p_sig2, p_p_bkg2}, {Asig2, Abkg2}, true); - - // Construct the subsidiary poisson measurements constraining the statistical fluctuations - RooArgSet globs2; - RooArgSet constraints2{createConstraints(gammas2, {p_h_sig, p_h_bkg}, globs2)}; - - // Construct the joint model - RooProdPdf model2("model2", "model2", {model2_tmp, constraints2}); - - // ************ Fit all models to data and plot ********************* - - std::unique_ptr result0{model0.fitTo(sumData, PrintLevel(0), Save())}; - std::unique_ptr result1{model1.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs1))}; - std::unique_ptr result2{model2.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs2))}; - - TCanvas *can = new TCanvas("can", "", 1500, 600); - can->Divide(3, 1); - - TPaveText pt(-19.5, 1, -2, 25); - pt.SetFillStyle(0); - pt.SetBorderSize(0); - - can->cd(1); - auto frame = x.frame(Title("No template uncertainties")); - // Plot data to enable automatic determination of model0 normalisation: - sumData.plotOn(frame); - model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0)); - // Plot data again to show it on top of model0 error bands: - sumData.plotOn(frame); - // Plot model components - model0.plotOn(frame, LineColor(kBlue)); - model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure)); - model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed)); - model0.paramOn(frame); - - sigData->plotOn(frame, MarkerColor(kBlue)); - frame->Draw(); - - for (auto text : {"No template uncertainties", "are taken into account.", "This leads to low errors", - "for the parameters A, since", "the only source of errors", "are the data statistics."}) { - pt.AddText(text); - } - pt.DrawClone(); - - can->cd(2); - frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately")); - sumData.plotOn(frame); - model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1)); - // Plot data again to show it on top of error bands: - sumData.plotOn(frame); - model1.plotOn(frame, LineColor(kBlue)); - model1.plotOn(frame, Components(p_p_sig1), LineColor(kAzure)); - model1.plotOn(frame, Components(p_p_bkg1), LineColor(kRed)); - model1.paramOn(frame, Parameters({Asig1, Abkg1})); - - sigData->plotOn(frame, MarkerColor(kBlue)); - frame->Draw(); - - pt.Clear(); - for (auto text : {"With gamma parameters, the", "signal & background templates", "can adapt to the data.", - "Note how the blue signal", "template changes its shape.", "This leads to higher errors", - "of the scale parameters A."}) { - pt.AddText(text); - } - pt.DrawClone(); - - can->cd(3); - frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)")); - sumData.plotOn(frame); - model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2)); - // Plot data again to show it on top of model0 error bands: - sumData.plotOn(frame); - model2.plotOn(frame, LineColor(kBlue)); - model2.plotOn(frame, Components(p_p_sig2), LineColor(kAzure)); - model2.plotOn(frame, Components(p_p_bkg2), LineColor(kRed)); - model2.paramOn(frame, Parameters({Asig2, Abkg2})); - - sigData->plotOn(frame, MarkerColor(kBlue)); - frame->Draw(); - - pt.Clear(); - for (auto text : {"When signal and background", "template share one gamma para-", "meter per bin, they adapt less.", - "The errors of the A parameters", "also shrink slightly."}) { - pt.AddText(text); - } - pt.DrawClone(); - - std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl; - std::cout << "Asig [BB ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl; - std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl; -} diff --git a/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py b/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py deleted file mode 100644 index a06d6876944cc..0000000000000 --- a/tutorials/roofit/roofit/rf709b_BarlowBeestonHistFactory.py +++ /dev/null @@ -1,251 +0,0 @@ -## \file -## \ingroup tutorial_roofit_main -## \notebook -js -## Implementing the Barlow-Beeston method for taking into account the statistical -## uncertainty of a Monte-Carlo fit template. -## -## \macro_image -## \macro_code -## \macro_output -## -## New implementation of the original rf709_BarlowBeeston tutorial, using the -## same kind of RooFit models that are built by HistFactory. Therefore, this -## tutorial is a good starting point if you want to understand how HistFactory -## works. -## -## \date November 2026 -## \author Jonas Rembser, CERN - -import math - -import ROOT - - -def create_constraints(gammas, samples, global_observables): - constraints = ROOT.RooArgSet() - for i, gamma in enumerate(gammas): - nom_name = "nom_" + gamma.GetName() - tau_name = gamma.GetName() + "_tau" - mean_name = gamma.GetName() + "_mean" - poisson_name = gamma.GetName() + "_constraint" - - nom_val = sum(func.dataHist().weight(i) for func in samples) - - # Like in HistFactory, the gamma of an empty bin is fixed and not - # constrained. Otherwise, it would be a flat direction of the likelihood. - # The other gammas get the relative Poisson uncertainty as initial error. - if nom_val == 0.0: - gamma.setConstant(True) - continue - gamma.setError(1.0 / math.sqrt(nom_val)) - - nominal = ROOT.RooRealVar(nom_name, nom_name, nom_val, 0.0, ROOT.RooNumber.infinity()) - nominal.setConstant(True) - global_observables.add(nominal) - - tau = ROOT.RooConstVar(tau_name, tau_name, nominal.getVal()) - - mean = ROOT.RooProduct(mean_name, mean_name, gamma, tau) - - poisson = ROOT.RooPoisson(poisson_name, poisson_name, nominal, mean) - - # The constraint pdf takes ownership of the objects it depends on, and - # the returned set takes ownership of the constraint pdf. - for helper in (nominal, tau, mean): - ROOT.SetOwnership(helper, False) - poisson.addOwnedComponents(ROOT.RooArgSet(nominal, tau, mean)) - - constraints.addOwned(poisson) - return constraints - - -ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.NumIntegration) - -# First, construct a likelihood model with a Gaussian signal on top of a uniform background -x = ROOT.RooRealVar("x", "x", -20, 20) -x.setBins(25) - -meanG = ROOT.RooRealVar("meanG", "meanG", 1, -10, 10) -sigG = ROOT.RooRealVar("sigG", "sigG", 1.5, -10, 10) -g = ROOT.RooGaussian("g", "Gauss", x, meanG, sigG) -u = ROOT.RooUniform("u", "Uniform", x) - -# Generate the data to be fitted -sigData = g.generate(x, 50) -bkgData = u.generate(x, 1000) - -sumData = ROOT.RooDataSet("sumData", "Gauss + Uniform", x) -sumData.append(sigData) -sumData.append(bkgData) - -# Make histogram templates for signal and background. -# Let's take a signal distribution with low statistics and a more accurate -# background distribution. -# Normally, these come from Monte Carlo simulations, but we will just generate them. -dh_sig = g.generateBinned(x, 50) -dh_bkg = u.generateBinned(x, 10000) - -# ***** Case 0 - 'Rigid templates' ***** - -# Construct histogram shapes for signal and background -p_h_sig = ROOT.RooHistFunc("p_h_sig", "p_h_sig", x, dh_sig) -p_h_bkg = ROOT.RooHistFunc("p_h_bkg", "p_h_bkg", x, dh_bkg) - -# Construct scale factors for adding the two distributions -Asig0 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) -Abkg0 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) - -# Construct the sum model -model0 = ROOT.RooRealSumPdf("model0", "model0", [p_h_sig, p_h_bkg], [Asig0, Abkg0], True) - -# ***** Case 1 - 'Barlow Beeston' ***** - -# Construct parameterized histogram shapes for signal and background -gammas_sig1 = ROOT.ParamHistFunc.createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10) -gammas_bkg1 = ROOT.ParamHistFunc.createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10) - -p_phf_sig1 = ROOT.ParamHistFunc("p_phf_sig1", "p_phf_sig1", [x], gammas_sig1) -p_phf_bkg1 = ROOT.ParamHistFunc("p_phf_bkg1", "p_phf_bkg1", [x], gammas_bkg1) - -p_p_sig1 = ROOT.RooProduct("p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1) -p_p_bkg1 = ROOT.RooProduct("p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1) - -Asig1 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) -Abkg1 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) - -# Construct the sum of these -model1_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig1, p_p_bkg1], [Asig1, Abkg1], True) - -# Construct the subsidiary poisson measurements constraining the histogram parameters -# These ensure that the bin contents of the histograms are only allowed to vary within -# the statistical uncertainty of the Monte Carlo. -globs1 = ROOT.RooArgSet() -constraints_sig1 = create_constraints(gammas_sig1, [p_h_sig], globs1) -constraints_bkg1 = create_constraints(gammas_bkg1, [p_h_bkg], globs1) - -# Construct the joint model with template PDFs and constraints -model1 = ROOT.RooProdPdf("model1", "model1", [model1_tmp, *constraints_sig1, *constraints_bkg1]) - -# ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) ***** - -# Construct the histogram shapes, using the same parameters for signal and background. -# There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to -# fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background -# histogram, so the bins fluctuate synchronously. This reduces the number of parameters. -gammas2 = ROOT.ParamHistFunc.createParamSet("p_ph_gamma", x.numBins(), 0, 10) - -p_phf2 = ROOT.ParamHistFunc("p_phf2", "p_phf2", [x], gammas2) - -p_p_sig2 = ROOT.RooProduct("p_p_sig2", "p_p_sig2", p_h_sig, p_phf2) -p_p_bkg2 = ROOT.RooProduct("p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2) - -Asig2 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) -Abkg2 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) - -# As before, construct the sum of signal2 and background2 -model2_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig2, p_p_bkg2], [Asig2, Abkg2], True) - -# Construct the subsidiary poisson measurements constraining the statistical fluctuations -globs2 = ROOT.RooArgSet() -constraints2 = create_constraints(gammas2, [p_h_sig, p_h_bkg], globs2) - -# Construct the joint model -model2 = ROOT.RooProdPdf("model2", "model2", [model2_tmp, *constraints2]) - -# ************ Fit all models to data and plot ********************* - -result0 = model0.fitTo(sumData, PrintLevel=0, Save=True) -result1 = model1.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs1) -result2 = model2.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs2) - -can = ROOT.TCanvas("can", "", 1500, 600) -can.Divide(3, 1) - -pt = ROOT.TPaveText(-19.5, 1, -2, 25) -pt.SetFillStyle(0) -pt.SetBorderSize(0) - -can.cd(1) -frame = x.frame(Title="No template uncertainties") -# Plot data to enable automatic determination of model0 normalisation: -sumData.plotOn(frame) -model0.plotOn(frame, LineColor="b", VisualizeError=result0) -# Plot data again to show it on top of model0 error bands: -sumData.plotOn(frame) -# Plot model components -model0.plotOn(frame, LineColor="b") -model0.plotOn(frame, Components=[p_h_sig], LineColor="kAzure") -model0.plotOn(frame, Components=[p_h_bkg], LineColor="r") -model0.paramOn(frame) - -sigData.plotOn(frame, MarkerColor="b") -frame.Draw() - -for text in [ - "No template uncertainties", - "are taken into account.", - "This leads to low errors", - "for the parameters A, since", - "the only source of errors", - "are the data statistics.", -]: - pt.AddText(text) -pt.DrawClone() - -can.cd(2) -frame = x.frame(Title="Barlow Beeston for Sig & Bkg separately") -sumData.plotOn(frame) -model1.plotOn(frame, LineColor="b", VisualizeError=result1) -# Plot data again to show it on top of error bands: -sumData.plotOn(frame) -model1.plotOn(frame, LineColor="b") -model1.plotOn(frame, Components=[p_p_sig1], LineColor="kAzure") -model1.plotOn(frame, Components=[p_p_bkg1], LineColor="r") -model1.paramOn(frame, Parameters=[Asig1, Abkg1]) - -sigData.plotOn(frame, MarkerColor="b") -frame.Draw() - -pt.Clear() -for text in [ - "With gamma parameters, the", - "signal & background templates", - "can adapt to the data.", - "Note how the blue signal", - "template changes its shape.", - "This leads to higher errors", - "of the scale parameters A.", -]: - pt.AddText(text) -pt.DrawClone() - -can.cd(3) -frame = x.frame(Title="Barlow Beeston light for (Sig+Bkg)") -sumData.plotOn(frame) -model2.plotOn(frame, LineColor="b", VisualizeError=result2) -# Plot data again to show it on top of model0 error bands: -sumData.plotOn(frame) -model2.plotOn(frame, LineColor="b") -model2.plotOn(frame, Components=[p_p_sig2], LineColor="kAzure") -model2.plotOn(frame, Components=[p_p_bkg2], LineColor="r") -model2.paramOn(frame, Parameters=[Asig2, Abkg2]) - -sigData.plotOn(frame, MarkerColor="b") -frame.Draw() - -pt.Clear() -for text in [ - "When signal and background", - "template share one gamma para-", - "meter per bin, they adapt less.", - "The errors of the A parameters", - "also shrink slightly.", -]: - pt.AddText(text) -pt.DrawClone() - -print(f"Asig [normal ] = {Asig0.getVal()} +/- {Asig0.getError()}") -print(f"Asig [BB ] = {Asig1.getVal()} +/- {Asig1.getError()}") -print(f"Asig [BBlight] = {Asig2.getVal()} +/- {Asig2.getError()}") - -can.SaveAs("rf709b_BarlowBeestonHistFactory.png") From e44e654fc264078467f88c140a303b5ff859ae0c Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 6 Nov 2026 12:39:30 +0000 Subject: [PATCH 5/5] [RF] Deprecate RooParamHistFunc and RooHistConstraint Their functionality is too close to the HistFactory classes, which has repeatedly caused confusion. Users should switch to ParamHistFunc with per-bin RooPoisson constraints, as shown in rf709_BarlowBeeston. Both classes will be removed in ROOT 6.46. The deprecation is not guarded by ROOFIT_BUILDS_ITSELF, because that definition also reaches CMake and interpreter users. Instead, the warnings are silenced locally where RooFit still uses the classes. Assisted-by: claude-code:claude-opus-5-5 --- README/ReleaseNotes/v644/index.md | 1 + roofit/codegen/inc/RooFit/CodegenImpl.h | 13 +++++++++++++ roofit/codegen/src/CodegenImpl.cxx | 13 +++++++++++++ roofit/roofit/inc/RooHistConstraint.h | 5 ++++- roofit/roofit/inc/RooParamHistFunc.h | 5 ++++- roofit/roofit/src/RooHistConstraint.cxx | 20 ++++++++++++++++++++ roofit/roofit/src/RooLagrangianMorphFunc.cxx | 13 +++++++++++++ roofit/roofit/src/RooParamHistFunc.cxx | 5 +++++ roofit/roofit/test/testRooParamHistFunc.cxx | 7 +++++++ 9 files changed, 80 insertions(+), 2 deletions(-) diff --git a/README/ReleaseNotes/v644/index.md b/README/ReleaseNotes/v644/index.md index 500a570b918a1..434d204570815 100644 --- a/README/ReleaseNotes/v644/index.md +++ b/README/ReleaseNotes/v644/index.md @@ -21,6 +21,7 @@ The following people have contributed to this new version: * The `RooTrace` class was untested and is now removed. For memory debugging, please use general-purpose tools such as AddressSanitizer or Valgrind instead. * Flattening nested **RooSimultaneous** in the constructor, deprecated in ROOT 6.42, is removed. RooSimultaneous components are now kept as they are, which also allows adding them with `RooSimultaneous::addPdf()`. They are only supported if their index category is a parameter ("switch mode"), which now also works for toy generation and for plotting slices of the outer RooSimultaneous. If the index category of a nested RooSimultaneous is in the dataset, creating the likelihood fails with a clear error for all evaluation backends. To combine channels, use one RooSimultaneous over a **RooCategory** with one state per channel, or over a **RooSuperCategory** if the channels span the full product of several categories. RooSimultaneous objects that were flattened by previous ROOT versions and are read from files are not affected. The internal `RooSimultaneous::InitializationOutput` struct and `RooSimultaneous::initialize()` function are removed as well. * The PyMVA library and the associated `tmva-pymva` build option that were deprecated in ROOT 6.42 have been removed +* The `RooParamHistFunc` and `RooHistConstraint` classes are deprecated and will be removed in ROOT 6.46. Their functionality is too close to the HistFactory classes, which has repeatedly caused confusion about which classes to use and how they differ. Please use the HistFactory class `ParamHistFunc` instead, multiplied with a `RooHistFunc` and constrained by one `RooPoisson` per bin, like the models built by HistFactory. The [rf709_BarlowBeeston](https://root.cern/doc/v644/rf709__BarlowBeeston_8C.html) tutorial was rewritten to show how this works. We are aware that these classes might have been used in RooWorkspaces that were written to disk. If you face problems because of that, please reach out to the RooFit developers by opening a [GitHub issue](https://github.com/root-project/root/issues) or posting on the [ROOT forum](https://root-forum.cern.ch/). * The **TMPIFile** class and the associated `mpi` build option (not to be confused with `minuit2_mpi`, which is unaffected) have been removed. diff --git a/roofit/codegen/inc/RooFit/CodegenImpl.h b/roofit/codegen/inc/RooFit/CodegenImpl.h index 9e02b392cf17a..eeccb859dc1db 100644 --- a/roofit/codegen/inc/RooFit/CodegenImpl.h +++ b/roofit/codegen/inc/RooFit/CodegenImpl.h @@ -116,7 +116,20 @@ void codegenImpl(RooLognormal &arg, CodegenContext &ctx); void codegenImpl(RooMultiPdf &arg, CodegenContext &ctx); void codegenImpl(RooMultiVarGaussian &arg, CodegenContext &ctx); void codegenImpl(RooONNXFunc &arg, CodegenContext &ctx); +// Support for the deprecated RooParamHistFunc, to be removed in ROOT 6.46 +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable : 4996) +#else +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wdeprecated-declarations" +#endif void codegenImpl(RooParamHistFunc &arg, CodegenContext &ctx); +#ifdef _MSC_VER +#pragma warning(pop) +#else +#pragma GCC diagnostic pop +#endif void codegenImpl(RooPoisson &arg, CodegenContext &ctx); void codegenImpl(RooPolyVar &arg, CodegenContext &ctx); void codegenImpl(RooPolynomial &arg, CodegenContext &ctx); diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index 3199a5311159a..e2b3fc3d99f8e 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -705,6 +705,14 @@ void codegenImpl(RooFit::Detail::RooNormalizedPdf &arg, CodegenContext &ctx) ctx.addResult(&arg, ctx.getResult(arg.pdf()) + "/" + ctx.getResult(arg.normIntegral())); } +// Support for the deprecated RooParamHistFunc, to be removed in ROOT 6.46 +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable : 4996) +#else +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wdeprecated-declarations" +#endif void codegenImpl(RooParamHistFunc &arg, CodegenContext &ctx) { std::string const &idx = arg.dataHist().calculateTreeIndexForCodeSquash(ctx, arg.xList()); @@ -720,6 +728,11 @@ void codegenImpl(RooParamHistFunc &arg, CodegenContext &ctx) } ctx.addResult(&arg, result.str()); } +#ifdef _MSC_VER +#pragma warning(pop) +#else +#pragma GCC diagnostic pop +#endif void codegenImpl(RooPoisson &arg, CodegenContext &ctx) { diff --git a/roofit/roofit/inc/RooHistConstraint.h b/roofit/roofit/inc/RooHistConstraint.h index 0e347cc940085..23a2b4c311b82 100644 --- a/roofit/roofit/inc/RooHistConstraint.h +++ b/roofit/roofit/inc/RooHistConstraint.h @@ -14,7 +14,10 @@ #include #include -class RooHistConstraint : public RooAbsPdf { +class R__DEPRECATED(6, 46, + "Its functionality is too close to the HistFactory classes. Please constrain the parameters of " + "the HistFactory class ParamHistFunc with one RooPoisson per bin instead, as shown in the " + "rf709_BarlowBeeston tutorial.") RooHistConstraint : public RooAbsPdf { public: RooHistConstraint() {} ; RooHistConstraint(const char *name, const char *title, const RooArgSet& phfSet, int threshold=1000000); diff --git a/roofit/roofit/inc/RooParamHistFunc.h b/roofit/roofit/inc/RooParamHistFunc.h index c227c89d70af1..ca1000eca4040 100644 --- a/roofit/roofit/inc/RooParamHistFunc.h +++ b/roofit/roofit/inc/RooParamHistFunc.h @@ -21,7 +21,10 @@ #include "RooDataHist.h" #include -class RooParamHistFunc : public RooAbsReal { +class R__DEPRECATED(6, 46, + "Its functionality is too close to the HistFactory classes. Please use the HistFactory class " + "ParamHistFunc instead, multiplying it with a RooHistFunc, as shown in the rf709_BarlowBeeston " + "tutorial.") RooParamHistFunc : public RooAbsReal { public: RooParamHistFunc() {} ; RooParamHistFunc(const char *name, const char *title, RooDataHist &dh, const RooAbsArg &x, diff --git a/roofit/roofit/src/RooHistConstraint.cxx b/roofit/roofit/src/RooHistConstraint.cxx index 039037383f141..4d8a76bd5c045 100644 --- a/roofit/roofit/src/RooHistConstraint.cxx +++ b/roofit/roofit/src/RooHistConstraint.cxx @@ -10,6 +10,11 @@ /** \class RooHistConstraint * \ingroup Roofit + * + * \deprecated RooHistConstraint will be removed in ROOT 6.46, because its functionality is too close to the + * HistFactory classes. Please constrain the parameters of the HistFactory class ParamHistFunc with one RooPoisson + * per bin instead, as shown in the rf709_BarlowBeeston tutorial. + * * The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties. * Following the Barlow-Beeston method, it adds Poisson constraints for each bin that * constrain the statistical uncertainty of the template histogram. @@ -25,6 +30,15 @@ #include +// RooHistConstraint is deprecated together with RooParamHistFunc, which it is +// built on, so its implementation necessarily refers to the deprecated class. +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable : 4996) +#else +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wdeprecated-declarations" +#endif //////////////////////////////////////////////////////////////////////////////// /// Create a new RooHistConstraint. @@ -217,3 +231,9 @@ double RooHistConstraint::getLogVal(const RooArgSet* /*set*/) const return sum ; } + +#ifdef _MSC_VER +#pragma warning(pop) +#else +#pragma GCC diagnostic pop +#endif diff --git a/roofit/roofit/src/RooLagrangianMorphFunc.cxx b/roofit/roofit/src/RooLagrangianMorphFunc.cxx index 4bdbc63796888..a4854630b0cf9 100644 --- a/roofit/roofit/src/RooLagrangianMorphFunc.cxx +++ b/roofit/roofit/src/RooLagrangianMorphFunc.cxx @@ -1610,6 +1610,14 @@ RooRealVar *RooLagrangianMorphFunc::setupObservable(const char *obsname, TClass obsExists = true; } +// Support for the deprecated RooParamHistFunc, to be removed in ROOT 6.46 +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable : 4996) +#else +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wdeprecated-declarations" +#endif if (mode && mode->InheritsFrom(RooHistFunc::Class())) { obs = static_cast(dynamic_cast(inputExample)->getHistObsList().first()); obsExists = true; @@ -1619,6 +1627,11 @@ RooRealVar *RooLagrangianMorphFunc::setupObservable(const char *obsname, TClass obsExists = true; _observables.add(*obs); } +#ifdef _MSC_VER +#pragma warning(pop) +#else +#pragma GCC diagnostic pop +#endif // Note: "found!" will be printed if s2 is a substring of s1, both s1 and s2 // are of type std::string. s1.find(s2) diff --git a/roofit/roofit/src/RooParamHistFunc.cxx b/roofit/roofit/src/RooParamHistFunc.cxx index 778ce475bdee3..27d7787c7ed78 100644 --- a/roofit/roofit/src/RooParamHistFunc.cxx +++ b/roofit/roofit/src/RooParamHistFunc.cxx @@ -10,6 +10,11 @@ /** \class RooParamHistFunc * \ingroup Roofit + * + * \deprecated RooParamHistFunc will be removed in ROOT 6.46, because its functionality is too close to the + * HistFactory classes. Please use the HistFactory class ParamHistFunc instead, multiplying it with a RooHistFunc, + * as shown in the rf709_BarlowBeeston tutorial. + * * A histogram function that assigns scale parameters to every bin. Instead of the bare bin contents, * it therefore yields: * \f[ diff --git a/roofit/roofit/test/testRooParamHistFunc.cxx b/roofit/roofit/test/testRooParamHistFunc.cxx index 262b978493f34..506ea26e9f720 100644 --- a/roofit/roofit/test/testRooParamHistFunc.cxx +++ b/roofit/roofit/test/testRooParamHistFunc.cxx @@ -14,6 +14,13 @@ #include +// RooParamHistFunc is deprecated, but still tested until its removal in ROOT 6.46 +#ifdef _MSC_VER +#pragma warning(disable : 4996) +#else +#pragma GCC diagnostic ignored "-Wdeprecated-declarations" +#endif + TEST(RooParamHistFunc, Integration) { // This tests the analytical integration of RooParamHistFunc,