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 80cfc3f7113d9..4d8a76bd5c045 100644 --- a/roofit/roofit/src/RooHistConstraint.cxx +++ b/roofit/roofit/src/RooHistConstraint.cxx @@ -10,14 +10,17 @@ /** \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. * * It can therefore be used to estimate the Monte Carlo uncertainty of a fit. * - * Check also the tutorial rf709_BarlowBeeston.C - * */ #include @@ -27,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. @@ -219,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 3f37cd548ece1..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[ @@ -19,8 +24,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/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, 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")