Skip to content

Commit

Permalink
[RF] Consider rangeName for integration in RooParamHistFunc
Browse files Browse the repository at this point in the history
The logic for summing over histogram bins in different ranges used in
RooHistPdf is also implemented in RooParamHistFunc. This means the
range is now considered when computing integrals of RooParamHistFunc.

RooParamHistFunc allows you to scale the counts in each bin with a
parameter. The interface of RooDataHist::sum was extended with a
function parameter to inject the logic of scaling the bin weight
depending on the bin index.

This commit partly fixes issue #7182. We still need to implement the
range feature in RooHistFunc.
  • Loading branch information
guitargeek committed Feb 19, 2021
1 parent b7e8bb7 commit c250ee4
Show file tree
Hide file tree
Showing 6 changed files with 61 additions and 28 deletions.
27 changes: 19 additions & 8 deletions roofit/roofit/src/RooParamHistFunc.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,11 @@
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include "RooRealVar.h"
#include "RooHelpers.h"
#include <math.h>
#include "TMath.h"


using namespace std ;

ClassImp(RooParamHistFunc);
Expand Down Expand Up @@ -244,19 +246,28 @@ Int_t RooParamHistFunc::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& a
/// Implement analytical integrations by doing appropriate weighting from component integrals
/// functions to integrators of components

Double_t RooParamHistFunc::analyticalIntegralWN(Int_t code, const RooArgSet* /*normSet2*/,const char* /*rangeName*/) const
Double_t RooParamHistFunc::analyticalIntegralWN(Int_t code, const RooArgSet* /*normSet2*/,const char* rangeName) const
{
// Supports only the scenario of integration over all dependents
R__ASSERT(code==1) ;

Double_t ret(0) ;
Int_t i(0) ;
for (const auto param : _p) {
auto p = static_cast<const RooAbsReal*>(param);
Double_t bin = p->getVal() ;
if (_relParam) bin *= getNominal(i++) ;
ret += bin ;
// The logic for summing over the histogram is borrowed from RooHistPdf with some differences:
//
// - a lambda function is used to inject the parameters for bin scaling into the RooDataHist::sum method
//
// - for simplicity, there is no check for the possibility of full-range integration with another overload of
// RooDataHist::sum
std::map<const RooAbsArg*, std::pair<Double_t, Double_t> > ranges;
for (const auto obs : _x) {
ranges[obs] = RooHelpers::getRangeOrBinningInterval(obs, rangeName);
}

auto getBinScale = [&](int iBin){ return static_cast<const RooAbsReal&>(_p[iBin]).getVal(); };

auto integrationSet = _dh.get();
RooArgSet sliceSet{};
Double_t ret = const_cast<RooDataHist&>(_dh).sum(*integrationSet, sliceSet, kFALSE, kTRUE, ranges, getBinScale);

// WVE fix this!!! Assume uniform binning for now!
Double_t binV(1) ;
for (const auto obs : _x) {
Expand Down
8 changes: 7 additions & 1 deletion roofit/roofitcore/inc/RooDataHist.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include <map>
#include <vector>
#include <string>
#include <functional>

class TObject ;
class RooAbsArg;
Expand Down Expand Up @@ -95,7 +96,12 @@ class RooDataHist : public RooAbsData, public RooDirItem {

Double_t sum(Bool_t correctForBinSize, Bool_t inverseCorr=kFALSE) const ;
Double_t sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, Bool_t correctForBinSize, Bool_t inverseCorr=kFALSE) ;
Double_t sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, Bool_t correctForBinSize, Bool_t inverseCorr, const std::map<const RooAbsArg*, std::pair<Double_t, Double_t> >& ranges);
Double_t sum(const RooArgSet& sumSet,
const RooArgSet& sliceSet,
Bool_t correctForBinSize,
Bool_t inverseCorr,
const std::map<const RooAbsArg*, std::pair<Double_t, Double_t> >& ranges,
std::function<double(int)> getBinScale = [](int){ return 1.0; } );

virtual Double_t weight() const {
// Return weight of current bin
Expand Down
4 changes: 4 additions & 0 deletions roofit/roofitcore/inc/RooHelpers.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@
#include <sstream>
#include <vector>
#include <string>
#include <utility>


namespace RooHelpers {

Expand Down Expand Up @@ -112,6 +114,8 @@ struct DisableCachingRAII {
};


std::pair<double, double> getRangeOrBinningInterval(RooAbsArg const* arg, const char* rangeName);


}

Expand Down
5 changes: 3 additions & 2 deletions roofit/roofitcore/src/RooDataHist.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -1553,7 +1553,8 @@ Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, Bo

Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet,
Bool_t correctForBinSize, Bool_t inverseBinCor,
const std::map<const RooAbsArg*, std::pair<Double_t, Double_t> >& ranges)
const std::map<const RooAbsArg*, std::pair<Double_t, Double_t> >& ranges,
std::function<double(int)> getBinScale)
{
checkInit();
checkBinBounds();
Expand Down Expand Up @@ -1624,7 +1625,7 @@ Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet,
const Double_t corr = correctForBinSize ? (inverseBinCor ? 1. / _binv[ibin] : _binv[ibin] ) : 1.0;
//cout << "adding bin[" << ibin << "] to sum wgt = " << _wgt[ibin] << " binv = " << theBinVolume << " _binv[" << ibin << "] " << _binv[ibin] << endl;
// const Double_t y = _wgt[ibin] * corr * corrPartial - carry;
const Double_t y = get_wgt(ibin) * corr * corrPartial - carry;
const Double_t y = getBinScale(ibin)*(get_wgt(ibin) * corr * corrPartial) - carry;
const Double_t t = total + y;
carry = (t - total) - y;
total = t;
Expand Down
26 changes: 26 additions & 0 deletions roofit/roofitcore/src/RooHelpers.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -144,4 +144,30 @@ void checkRangeOfParameters(const RooAbsReal* callingClass, std::initializer_lis
}
}


/// Get the lower and upper bound of parameter range if arg can be casted to RooAbsRealLValue.
/// If no range with rangeName is defined for the argument, this will check if a binning of the
/// same name exists and return the interval covered by the binning.
/// Returns `{-infinity, infinity}` if agument can't be casted to RooAbsRealLValue* or if no
/// range or binning with the requested name exists.
/// \param[in] arg RooAbsArg for which to get the range.
/// \param[in] rangeName The name of the range.
std::pair<double, double> getRangeOrBinningInterval(RooAbsArg const* arg, const char* rangeName) {
auto rlv = dynamic_cast<RooAbsRealLValue const*>(arg);
if (rlv) {
const RooAbsBinning* binning = rlv->getBinningPtr(rangeName);
if (rangeName && rlv->hasRange(rangeName)) {
return {rlv->getMin(rangeName), rlv->getMax(rangeName)};
} else if (binning) {
if (!binning->isParameterized()) {
return {binning->lowBound(), binning->highBound()};
} else {
return {binning->lowBoundFunc()->getVal(), binning->highBoundFunc()->getVal()};
}
}
}
return {-std::numeric_limits<double>::infinity(), +std::numeric_limits<double>::infinity()};
}


}
19 changes: 2 additions & 17 deletions roofit/roofitcore/src/RooHistPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ discrete dimensions.
#include "RooCategory.h"
#include "RooWorkspace.h"
#include "RooGlobalFunc.h"
#include "RooHelpers.h"

#include "TError.h"
#include "TBuffer.h"
Expand Down Expand Up @@ -317,7 +318,6 @@ Int_t RooHistPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars,
}



////////////////////////////////////////////////////////////////////////////////
/// Return integral identified by 'code'. The actual integration
/// is deferred to RooDataHist::sum() which implements partial
Expand All @@ -342,22 +342,7 @@ Double_t RooHistPdf::analyticalIntegral(Int_t code, const char* rangeName) const
intSet.add(*ha);
}
if (!(code & 1)) {
RooAbsRealLValue* rlv = dynamic_cast<RooAbsRealLValue*>(pa);
if (rlv) {
const RooAbsBinning* binning = rlv->getBinningPtr(rangeName);
if (rangeName && rlv->hasRange(rangeName)) {
ranges[ha] = std::make_pair(
rlv->getMin(rangeName), rlv->getMax(rangeName));
} else if (binning) {
if (!binning->isParameterized()) {
ranges[ha] = std::make_pair(
binning->lowBound(), binning->highBound());
} else {
ranges[ha] = std::make_pair(
binning->lowBoundFunc()->getVal(), binning->highBoundFunc()->getVal());
}
}
}
ranges[ha] = RooHelpers::getRangeOrBinningInterval(pa, rangeName);
}
// WVE must sync hist slice list values to pdf slice list
// Transfer values from
Expand Down

1 comment on commit c250ee4

@TheQuantiser
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hello, I find that adding "range" doesn't change the integral error. I have posted my issue on root forums [*]. Could you please take a look?

[*] https://root-forum.cern.ch/t/barlow-beeston-in-subrange/43909

Please sign in to comment.