From 5f027832bb708e1735dd19c3d962d2e4a6b82981 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sat, 13 Apr 2024 02:56:30 +0200 Subject: [PATCH] [RF] Fix computation graphs with RooParamHistFunc The `RooParamHistFunc` doesn't take any observable RooRealVar as constructor argument. It assumes as observable the internal variables in the passed RooDataHist. This means it is in most contexts unusable, because the input can't be changed, other than loading a different bin in the dataset. This also breaks the Barlow-Beeston tutorial since the new evaluation backend is the default, and it is more sentitive to these issues: https://root.cern.ch/doc/master/rf709__BarlowBeeston_8C.html There was actually a constructor that took a `RooAbsArg x`, but it was simply ignored. To fix all these problems, the existing constructors were replaced by a new one that takes the observable explicitly, and this is mentioned in the release notes. The class is not used much, because ussually people use HistFactory or CMS combine for these kind of fits. --- README/ReleaseNotes/v632/index.md | 17 +++ roofit/roofit/inc/RooParamHistFunc.h | 5 +- roofit/roofit/src/RooParamHistFunc.cxx | 114 ++++++-------------- roofit/roofit/test/testRooParamHistFunc.cxx | 4 +- tutorials/roofit/rf709_BarlowBeeston.C | 8 +- tutorials/roofit/rf709_BarlowBeeston.py | 8 +- 6 files changed, 59 insertions(+), 97 deletions(-) diff --git a/README/ReleaseNotes/v632/index.md b/README/ReleaseNotes/v632/index.md index b20e8df7914e4..fd8fb172f93af 100644 --- a/README/ReleaseNotes/v632/index.md +++ b/README/ReleaseNotes/v632/index.md @@ -223,6 +223,23 @@ Instantiating the following classes and even including their header files is dep Please use the higher-level functions `RooAbsPdf::createNLL()` and `RooAbsPdf::createChi2()` if you want to create objects that represent test statistics. +### Change of RooParamHistFunc + +The `RooParamHistFunc` didn't take any observable `RooRealVar` as constructor +argument. It assumes as observable the internal variables in the passed +RooDataHist. This means it was in most contexts unusable, because the input +can't be changed, other than loading a different bin in the dataset. + +Furthermore, there was actually a constructor that took a `RooAbsArg x`, but it +was simply ignored. + +To fix all these problems, the existing constructors were replaced by a new one +that takes the observable explicitly. + +Since the old constructors resulted in wrong computation graphs that caused +trouble with the new CPU evaluation backend, they had to be removed without +deprecation. Please adapt your code if necessary. + ## RDataFrame * The RDataFrame constructors that take in input one or more file names (or globs thereof) will now infer the format of the dataset, either TTree or RNTuple, that is stored in the first input file. When multiple files are specified, it is assumed that all other files contain a coherent dataset of the same format and with the same schema, exactly as it used to happen with TChain. This automatic inference further contributes towards a zero-code-change experience when moving from processing a TTree to processing an RNTuple dataset while using an RDataFrame. It also introduces a backwards-incompatible behaviour, i.e. now the constructor needs to open one file in order to infer the dataset type. This means that if the file does not exist, the constructor will throw an exception. Previously, an exception would be thrown only at a JIT-ting time, before the start of the computations. diff --git a/roofit/roofit/inc/RooParamHistFunc.h b/roofit/roofit/inc/RooParamHistFunc.h index 97e07c95fdaeb..b42392bd59e33 100644 --- a/roofit/roofit/inc/RooParamHistFunc.h +++ b/roofit/roofit/inc/RooParamHistFunc.h @@ -24,9 +24,8 @@ class RooParamHistFunc : public RooAbsReal { public: RooParamHistFunc() {} ; - RooParamHistFunc(const char *name, const char *title, RooDataHist& dh, bool paramRelative=true); - RooParamHistFunc(const char *name, const char *title, const RooAbsArg& x, RooDataHist& dh, bool paramRelative=true); - RooParamHistFunc(const char *name, const char *title, RooDataHist& dh, const RooParamHistFunc& paramSource, bool paramRelative=true) ; + RooParamHistFunc(const char *name, const char *title, RooDataHist &dh, const RooAbsArg &x, + const RooParamHistFunc *paramSource = nullptr, bool paramRelative = true); RooParamHistFunc(const RooParamHistFunc& other, const char* name=nullptr) ; TObject* clone(const char* newname) const override { return new RooParamHistFunc(*this,newname); } diff --git a/roofit/roofit/src/RooParamHistFunc.cxx b/roofit/roofit/src/RooParamHistFunc.cxx index 187e48b76201f..d637c53e416b6 100644 --- a/roofit/roofit/src/RooParamHistFunc.cxx +++ b/roofit/roofit/src/RooParamHistFunc.cxx @@ -23,88 +23,38 @@ * See also the tutorial rf709_BarlowBeeston.C */ -#include "Riostream.h" -#include "RooParamHistFunc.h" -#include "RooAbsCategory.h" -#include "RooRealVar.h" -#include "RooFitImplHelpers.h" -#include -#include "TMath.h" - - -using std::list; +#include +#include +#include ClassImp(RooParamHistFunc); -//////////////////////////////////////////////////////////////////////////////// -/// Populate x with observables - -RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataHist& dh, bool paramRelative) : - RooAbsReal(name,title), - _x("x","x",this), - _p("p","p",this), - _dh(dh), - _relParam(paramRelative) -{ - _x.add(*_dh.get()) ; - - // Now populate p with parameters - RooArgSet allVars ; - for (Int_t i=0 ; i<_dh.numEntries() ; i++) { - _dh.get(i) ; - - const char* vname = Form("%s_gamma_bin_%i",GetName(),i) ; - RooRealVar* var = new RooRealVar(vname,vname,0,1000) ; - var->setVal(_relParam ? 1 : _dh.weight()) ; - var->setError(_relParam ? 1 / sqrt(_dh.weight()) : sqrt(_dh.weight())); - var->setConstant(true) ; - allVars.add(*var) ; - _p.add(*var) ; - } - addOwnedComponents(allVars) ; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Populate x with observables - -RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, const RooAbsArg& /*x*/, RooDataHist& dh, bool paramRelative) : - RooAbsReal(name,title), - _x("x","x",this), - _p("p","p",this), - _dh(dh), - _relParam(paramRelative) +RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataHist &dh, const RooAbsArg &x, + const RooParamHistFunc *paramSource, bool paramRelative) + : RooAbsReal(name, title), _x("x", "x", this), _p("p", "p", this), _dh(dh), _relParam(paramRelative) { - _x.add(*_dh.get()) ; - - // Now populate p with parameters - RooArgSet allVars ; - for (Int_t i=0 ; i<_dh.numEntries() ; i++) { - _dh.get(i) ; - const char* vname = Form("%s_gamma_bin_%i",GetName(),i) ; - RooRealVar* var = new RooRealVar(vname,vname,0,1000) ; - var->setVal(_relParam ? 1 : _dh.weight()) ; - var->setError(_relParam ? 1 / sqrt(_dh.weight()) : sqrt(_dh.weight())); - var->setConstant(true) ; - allVars.add(*var) ; - _p.add(*var) ; - } - addOwnedComponents(allVars) ; -} - -//////////////////////////////////////////////////////////////////////////////// + // Populate x with observables + _x.add(x); -RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataHist& dh, const RooParamHistFunc& paramSource, bool paramRelative) : - RooAbsReal(name,title), - _x("x","x",this), - _p("p","p",this), - _dh(dh), - _relParam(paramRelative) -{ - // Populate x with observables - _x.add(*_dh.get()) ; + if (paramSource) { + // Now populate p with existing parameters + _p.add(paramSource->_p); + return; + } - // Now populate p with existing parameters - _p.add(paramSource._p) ; + // Now populate p with parameters + RooArgSet allVars; + for (Int_t i = 0; i < _dh.numEntries(); i++) { + _dh.get(i); + const char *vname = Form("%s_gamma_bin_%i", GetName(), i); + RooRealVar *var = new RooRealVar(vname, vname, 0, 1000); + var->setVal(_relParam ? 1 : _dh.weight()); + var->setError(_relParam ? 1 / sqrt(_dh.weight()) : sqrt(_dh.weight())); + var->setConstant(true); + allVars.add(*var); + _p.add(*var); + } + addOwnedComponents(allVars); } //////////////////////////////////////////////////////////////////////////////// @@ -124,11 +74,7 @@ double RooParamHistFunc::evaluate() const { Int_t idx = ((RooDataHist&)_dh).getIndex(_x,true) ; double ret = (static_cast(_p.at(idx)))->getVal() ; - if (_relParam) { - double nom = getNominal(idx) ; - ret *= nom ; - } - return ret ; + return _relParam ? ret * getNominal(idx) : ret; } void RooParamHistFunc::translate(RooFit::Detail::CodeSquashContext &ctx) const @@ -185,7 +131,7 @@ double RooParamHistFunc::getNominalError(Int_t ibin) const /// as the recursive division strategy of RooCurve cannot deal efficiently /// with the vertical lines that occur in a non-interpolated histogram -list* RooParamHistFunc::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const +std::list* RooParamHistFunc::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const { // Check that observable is in dataset, if not no hint is generated RooAbsLValue* lvarg = dynamic_cast(_dh.get()->find(obs.GetName())) ; @@ -197,7 +143,7 @@ list* RooParamHistFunc::plotSamplingHint(RooAbsRealLValue& obs, double x const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr); double* boundaries = binning->array() ; - list* hint = new list ; + std::list* hint = new std::list ; // Widen range slightly xlo = xlo - 0.01*(xhi-xlo) ; @@ -234,7 +180,7 @@ std::list* RooParamHistFunc::binBoundaries(RooAbsRealLValue& obs, double const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr); double* boundaries = binning->array() ; - list* hint = new list ; + std::list* hint = new std::list ; // Construct array with pairs of points positioned epsilon to the left and // right of the bin boundaries diff --git a/roofit/roofit/test/testRooParamHistFunc.cxx b/roofit/roofit/test/testRooParamHistFunc.cxx index 5244ceb6f201e..262b978493f34 100644 --- a/roofit/roofit/test/testRooParamHistFunc.cxx +++ b/roofit/roofit/test/testRooParamHistFunc.cxx @@ -35,7 +35,7 @@ TEST(RooParamHistFunc, Integration) } RooDataHist dh("dh", "dh", x, &h1); - RooParamHistFunc phf("phf", "", x, dh); + RooParamHistFunc phf("phf", "", dh, x); x.setRange("R1", 0, xMax * 0.5); std::unique_ptr integral{phf.createIntegral(x, x)}; @@ -87,7 +87,7 @@ TEST(RooParamHistFunc, IntegrationAndCloning) h1.FillRandom("f1", 50); RooDataHist dh1("dh1", "dh1", x, &h1); - RooParamHistFunc ph("ph", "", x, dh1); + RooParamHistFunc ph("ph", "", dh1, x); // Combine the RooParamHistFunc with something else in a RooRealSumPdf. // This is do make the test more similar to the Barlow-Beeston test, diff --git a/tutorials/roofit/rf709_BarlowBeeston.C b/tutorials/roofit/rf709_BarlowBeeston.C index fd8c2a78297c3..3637ffb5658dd 100644 --- a/tutorials/roofit/rf709_BarlowBeeston.C +++ b/tutorials/roofit/rf709_BarlowBeeston.C @@ -83,8 +83,8 @@ void rf709_BarlowBeeston() // ***** Case 1 - 'Barlow Beeston' ***** // Construct parameterized histogram shapes for signal and background - RooParamHistFunc p_ph_sig1("p_ph_sig","p_ph_sig",*dh_sig); - RooParamHistFunc p_ph_bkg1("p_ph_bkg","p_ph_bkg",*dh_bkg); + 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); RooRealVar Asig1("Asig","Asig",1,0.01,5000); RooRealVar Abkg1("Abkg","Abkg",1,0.01,5000); @@ -115,8 +115,8 @@ void rf709_BarlowBeeston() // 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); - RooParamHistFunc p_ph_bkg2("p_ph_bkg2", "p_ph_bkg2", *dh_bkg, p_ph_sig2, true); + 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); diff --git a/tutorials/roofit/rf709_BarlowBeeston.py b/tutorials/roofit/rf709_BarlowBeeston.py index 0becaec859e21..f0373a9a8577a 100644 --- a/tutorials/roofit/rf709_BarlowBeeston.py +++ b/tutorials/roofit/rf709_BarlowBeeston.py @@ -59,8 +59,8 @@ # 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) -p_ph_bkg1 = ROOT.RooParamHistFunc("p_ph_bkg", "p_ph_bkg", dh_bkg) +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) Asig1 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) Abkg1 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000) @@ -87,8 +87,8 @@ # 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) -p_ph_bkg2 = ROOT.RooParamHistFunc("p_ph_bkg2", "p_ph_bkg2", dh_bkg, p_ph_sig2, True) +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) Asig2 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000) Abkg2 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000)