diff --git a/roofit/roofit/inc/RooParametricStepFunction.h b/roofit/roofit/inc/RooParametricStepFunction.h index a375e2d5a9118..48defbbc1e627 100644 --- a/roofit/roofit/inc/RooParametricStepFunction.h +++ b/roofit/roofit/inc/RooParametricStepFunction.h @@ -15,12 +15,12 @@ #ifndef ROO_PARAMETRIC_STEP_FUNCTION #define ROO_PARAMETRIC_STEP_FUNCTION -#include "TArrayD.h" -#include "RooAbsPdf.h" -#include "RooRealProxy.h" -#include "RooListProxy.h" +#include +#include +#include + +#include -class RooRealVar; class RooArgList ; class RooParametricStepFunction : public RooAbsPdf { @@ -29,7 +29,7 @@ class RooParametricStepFunction : public RooAbsPdf { RooParametricStepFunction() {} RooParametricStepFunction(const char *name, const char *title, - RooAbsReal& x, const RooArgList& coefList, TArrayD& limits, Int_t nBins=1) ; + RooAbsReal& x, const RooArgList& coefList, TArrayD const& limits, Int_t nBins=1) ; RooParametricStepFunction(const RooParametricStepFunction& other, const char* name = nullptr); TObject* clone(const char* newname) const override { return new RooParametricStepFunction(*this, newname); } @@ -39,6 +39,8 @@ class RooParametricStepFunction : public RooAbsPdf { Int_t getnBins() const { return _nBins; } double* getLimits() { return _limits.GetArray(); } + std::list* plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const override ; + protected: double lastBinValue() const ; diff --git a/roofit/roofit/inc/RooStepFunction.h b/roofit/roofit/inc/RooStepFunction.h index 6c0d7f08a8d13..879e660a88904 100644 --- a/roofit/roofit/inc/RooStepFunction.h +++ b/roofit/roofit/inc/RooStepFunction.h @@ -17,18 +17,16 @@ #ifndef ROO_STEP_FUNCTION #define ROO_STEP_FUNCTION -#include "TArrayD.h" -#include "RooAbsReal.h" -#include "RooRealProxy.h" -#include "RooListProxy.h" +#include +#include +#include -class RooRealVar; class RooArgList ; class RooStepFunction : public RooAbsReal { public: - RooStepFunction() ; + RooStepFunction() {} RooStepFunction(const char *name, const char *title, RooAbsReal& x, const RooArgList& coefList, const RooArgList& limits, bool interpolate=false) ; @@ -38,6 +36,8 @@ class RooStepFunction : public RooAbsReal { const RooArgList& coefficients() { return _coefList; } const RooArgList& boundaries() { return _boundaryList; } + std::list* plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const override ; + protected: double evaluate() const override; @@ -47,7 +47,7 @@ class RooStepFunction : public RooAbsReal { RooRealProxy _x; RooListProxy _coefList ; RooListProxy _boundaryList ; - bool _interpolate ; + bool _interpolate = false; ClassDefOverride(RooStepFunction,1) // Step Function }; diff --git a/roofit/roofit/src/RooParametricStepFunction.cxx b/roofit/roofit/src/RooParametricStepFunction.cxx index 7579f53d31263..229983ae1fc81 100644 --- a/roofit/roofit/src/RooParametricStepFunction.cxx +++ b/roofit/roofit/src/RooParametricStepFunction.cxx @@ -78,15 +78,10 @@ xframe->Draw(); ~~~ */ -#include "Riostream.h" -#include "TArrayD.h" -#include +#include -#include "RooParametricStepFunction.h" -#include "RooRealVar.h" -#include "RooArgList.h" - -#include "TError.h" +#include +#include ClassImp(RooParametricStepFunction); @@ -94,7 +89,7 @@ ClassImp(RooParametricStepFunction); /// Constructor RooParametricStepFunction::RooParametricStepFunction(const char* name, const char* title, - RooAbsReal& x, const RooArgList& coefList, TArrayD& limits, Int_t nBins) : + RooAbsReal& x, const RooArgList& coefList, TArrayD const& limits, Int_t nBins) : RooAbsPdf(name, title), _x("x", "Dependent", this, x), _coefList("coefList","List of coefficients",this), @@ -110,9 +105,11 @@ RooParametricStepFunction::RooParametricStepFunction(const char* name, const cha for (auto *coef : coefList) { if (!dynamic_cast(coef)) { - std::cout << "RooParametricStepFunction::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName() - << " is not of type RooAbsReal" << std::endl ; - R__ASSERT(0) ; + std::stringstream errorMsg; + errorMsg << "RooParametricStepFunction::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName() + << " is not of type RooAbsReal"; + coutE(InputArguments) << errorMsg.str() << std::endl; + throw std::invalid_argument(errorMsg.str().c_str()); } _coefList.add(*coef) ; } @@ -144,10 +141,8 @@ Int_t RooParametricStepFunction::getAnalyticalIntegral(RooArgSet& allVars, RooAr //////////////////////////////////////////////////////////////////////////////// -double RooParametricStepFunction::analyticalIntegral(Int_t code, const char* rangeName) const +double RooParametricStepFunction::analyticalIntegral(Int_t /*code*/, const char* rangeName) const { - R__ASSERT(code==1) ; - // Case without range is trivial: p.d.f is by construction normalized if (!rangeName) { return 1.0 ; @@ -210,20 +205,10 @@ double RooParametricStepFunction::evaluate() const // in Bin i-1 (starting with Bin 0) if (i<_nBins) { // not in last Bin - RooRealVar* tmp = (RooRealVar*) _coefList.at(i-1); - value = tmp->getVal(); + value = static_cast(_coefList.at(i-1))->getVal(); break; } else { - // in last Bin - double sum(0.); - double binSize(0.); - for (Int_t j=1;j<_nBins;j++){ - RooRealVar* tmp = (RooRealVar*) _coefList.at(j-1); - binSize = _limits[j] - _limits[j-1]; - sum = sum + tmp->getVal()*binSize; - } - binSize = _limits[_nBins] - _limits[_nBins-1]; - value = (1.0 - sum)/binSize; + value = lastBinValue(); if (value<=0.0){ value = 0.000001; // cout << "RooParametricStepFunction: sum of values gt 1.0 -- beware!!" <* RooParametricStepFunction::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const +{ + if(obs.namePtr() != _x->namePtr()) { + return nullptr; + } + + // Retrieve position of all bin boundaries + std::span boundaries{_limits.GetArray(), static_cast(_limits.GetSize())}; + + // Use the helper function from RooCurve to make sure to get sampling hints + // that work with the RooFitPlotting. + return RooCurve::plotSamplingHintForBinBoundaries(boundaries, xlo, xhi); +} diff --git a/roofit/roofit/src/RooStepFunction.cxx b/roofit/roofit/src/RooStepFunction.cxx index d3be67f2e0eb2..b70b20a6c00b4 100644 --- a/roofit/roofit/src/RooStepFunction.cxx +++ b/roofit/roofit/src/RooStepFunction.cxx @@ -28,31 +28,19 @@ Note that in contrast to RooParametricStepFunction, a RooStepFunction is NOT a P but a not-normalized function (RooAbsReal) **/ -#include "Riostream.h" -#include "TArrayD.h" -#include +#include -#include "RooStepFunction.h" -#include "RooRealVar.h" -#include "RooArgList.h" -#include "RooMsgService.h" -#include "RooMath.h" - -using namespace std; +#include +#include +#include +#include +#include ClassImp(RooStepFunction); //////////////////////////////////////////////////////////////////////////////// /// Constructor -RooStepFunction::RooStepFunction() -{ - _interpolate = false ; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Constructor - RooStepFunction::RooStepFunction(const char* name, const char* title, RooAbsReal& x, const RooArgList& coefList, const RooArgList& boundaryList, bool interpolate) : RooAbsReal(name, title), @@ -63,8 +51,8 @@ RooStepFunction::RooStepFunction(const char* name, const char* title, { for (auto *coef : coefList) { if (!dynamic_cast(coef)) { - cout << "RooStepFunction::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName() - << " is not of type RooAbsReal" << endl ; + std::cout << "RooStepFunction::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName() + << " is not of type RooAbsReal" << std::endl ; assert(0) ; } _coefList.add(*coef) ; @@ -72,16 +60,16 @@ RooStepFunction::RooStepFunction(const char* name, const char* title, for (auto *boundary : boundaryList) { if (!dynamic_cast(boundary)) { - cout << "RooStepFunction::ctor(" << GetName() << ") ERROR: boundary " << boundary->GetName() - << " is not of type RooAbsReal" << endl ; + std::cout << "RooStepFunction::ctor(" << GetName() << ") ERROR: boundary " << boundary->GetName() + << " is not of type RooAbsReal" << std::endl; assert(0) ; } _boundaryList.add(*boundary) ; } - if (_boundaryList.getSize()!=_coefList.getSize()+1) { - coutE(InputArguments) << "RooStepFunction::ctor(" << GetName() << ") ERROR: Number of boundaries must be number of coefficients plus 1" << endl ; - throw string("RooStepFunction::ctor() ERROR: Number of boundaries must be number of coefficients plus 1") ; + if (_boundaryList.size()!=_coefList.size()+1) { + coutE(InputArguments) << "RooStepFunction::ctor(" << GetName() << ") ERROR: Number of boundaries must be number of coefficients plus 1" << std::endl ; + throw std::invalid_argument("RooStepFunction::ctor() ERROR: Number of boundaries must be number of coefficients plus 1") ; } } @@ -100,12 +88,12 @@ RooStepFunction::RooStepFunction(const RooStepFunction& other, const char* name) //////////////////////////////////////////////////////////////////////////////// -/// Transfer contents to vector for use below +/// Transfer contents to std::vector for use below double RooStepFunction::evaluate() const { - vector b(_boundaryList.getSize()) ; - vector c(_coefList.getSize()+3) ; + std::vector b(_boundaryList.size()) ; + std::vector c(_coefList.size()+3) ; Int_t nb(0) ; for (auto * boundary : static_range_cast(_boundaryList)) { b[nb++] = boundary->getVal() ; @@ -136,7 +124,7 @@ double RooStepFunction::evaluate() const // Make array of (0,coefficient values,0) Int_t nc(0) ; - vector y(_coefList.size()+3) ; + std::vector y(_coefList.size()+3) ; y[nc++] = 0 ; for(auto * coef : static_range_cast(_coefList)) { y[nc++] = coef->getVal() ; @@ -153,3 +141,22 @@ double RooStepFunction::evaluate() const return 0; } } + + +std::list *RooStepFunction::plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const +{ + if (obs.namePtr() != _x->namePtr()) { + return nullptr; + } + + // Retrieve position of all bin boundaries + std::vector boundaries; + boundaries.reserve(_boundaryList.size()); + for (auto *boundary : static_range_cast(_boundaryList)) { + boundaries.push_back(boundary->getVal()); + } + + // Use the helper function from RooCurve to make sure to get sampling hints + // that work with the RooFitPlotting. + return RooCurve::plotSamplingHintForBinBoundaries(boundaries, xlo, xhi); +} diff --git a/roofit/roofitcore/inc/RooCurve.h b/roofit/roofitcore/inc/RooCurve.h index 043f96721c307..e2c657b63ccd5 100644 --- a/roofit/roofitcore/inc/RooCurve.h +++ b/roofit/roofitcore/inc/RooCurve.h @@ -16,11 +16,15 @@ #ifndef ROO_CURVE #define ROO_CURVE -#include "TGraph.h" -#include "RooPlotable.h" +#include + +#include + +#include +#include + #include #include -#include "TMatrixDfwd.h" class RooAbsReal; class RooRealVar; @@ -71,6 +75,8 @@ class RooCurve : public TGraph, public RooPlotable { RooCurve* makeErrorBand(const std::vector& variations, double Z=1) const ; RooCurve* makeErrorBand(const std::vector& plusVar, const std::vector& minusVar, const TMatrixD& V, double Z=1) const ; + static std::list* plotSamplingHintForBinBoundaries(std::span boundaries, double xlo, double xhi); + protected: void calcBandInterval(const std::vector& variations,Int_t i,double Z,double& lo, double& hi, bool approxGauss) const ; @@ -83,13 +89,20 @@ class RooCurve : public TGraph, public RooPlotable { Int_t numee=0, bool doEEVal=false, double eeVal=0.0,std::list* samplingHint=nullptr) ; void addRange(const RooAbsFunc& func, double x1, double x2, double y1, double y2, double minDy, double minDx, - Int_t numee=0, bool doEEVal=false, double eeVal=0.0); + int numee, bool doEEVal, double eeVal, double epsilon); void shiftCurveToZero(); bool _showProgress = false; /// yval(minPoints); // Get list of initial x values. If function provides sampling hint use that, @@ -369,7 +370,7 @@ void RooCurve::addPoints(const RooAbsFunc &func, double xlo, double xhi, // If precision is <0, no attempt at recursive interpolation is made addPoint(x2,yval[step]) ; } else { - addRange(func,x1,x2,yval[step-1],yval[step],prec*yrangeEst,minDx,numee,doEEVal,eeVal); + addRange(func,x1,x2,yval[step-1],yval[step],prec*yrangeEst,minDx,numee,doEEVal,eeVal,epsilon); } step++ ; } @@ -394,10 +395,10 @@ void RooCurve::addPoints(const RooAbsFunc &func, double xlo, double xhi, void RooCurve::addRange(const RooAbsFunc& func, double x1, double x2, double y1, double y2, double minDy, double minDx, - Int_t numee, bool doEEVal, double eeVal) + int numee, bool doEEVal, double eeVal, double epsilon) { // Explicitly skip empty ranges to eliminate point duplication - if (std::abs(x2-x1)<1e-20) { + if (std::abs(x2-x1) <= epsilon) { return ; } @@ -424,8 +425,8 @@ void RooCurve::addRange(const RooAbsFunc& func, double x1, double x2, double dy= ymid - 0.5*(y1+y2); if((xmid - x1 >= minDx) && std::abs(dy)>0 && std::abs(dy) >= minDy) { // fill in each subrange - addRange(func,x1,xmid,y1,ymid,minDy,minDx,numee,doEEVal,eeVal); - addRange(func,xmid,x2,ymid,y2,minDy,minDx,numee,doEEVal,eeVal); + addRange(func,x1,xmid,y1,ymid,minDy,minDx,numee,doEEVal,eeVal,epsilon); + addRange(func,xmid,x2,ymid,y2,minDy,minDx,numee,doEEVal,eeVal,epsilon); } else { // add the endpoint @@ -886,3 +887,37 @@ bool RooCurve::isIdentical(const RooCurve& other, double tol, bool verbose) cons } + +//////////////////////////////////////////////////////////////////////////////// +/// Returns sampling hints for a histogram with given boundaries. This helper +/// function is meant to be used by binned RooAbsReals to produce sampling +/// hints that are working well with RooFits plotting. + +std::list * +RooCurve::plotSamplingHintForBinBoundaries(std::span boundaries, double xlo, double xhi) +{ + auto hint = new std::list; + + // Make sure the difference between two points around a bin boundary is + // larger than the relative epsilon for which the RooCurve considers two + // points as the same. Otherwise, the points right of the bin boundary would + // be skipped. + const double delta = (xhi - xlo) * RooCurve::relativeXEpsilon(); + + // Sample points right next to the plot limits + hint->push_back(xlo + delta); + hint->push_back(xhi - delta); + + // Sample points very close to the left and right of the bin boundaries that + // are strictly in between the plot limits. + for (const double x : boundaries) { + if (x - xlo > delta && xhi - x > delta) { + hint->push_back(x - delta); + hint->push_back(x + delta); + } + } + + hint->sort(); + + return hint; +} diff --git a/roofit/roofitcore/src/RooHistPdf.cxx b/roofit/roofitcore/src/RooHistPdf.cxx index 1a742ae244fbf..9a63e01eb63c5 100644 --- a/roofit/roofitcore/src/RooHistPdf.cxx +++ b/roofit/roofitcore/src/RooHistPdf.cxx @@ -28,6 +28,7 @@ discrete dimensions. #include "Riostream.h" #include "RooHistPdf.h" +#include "RooCurve.h" #include "RooDataHist.h" #include "RooMsgService.h" #include "RooRealVar.h" @@ -511,26 +512,11 @@ std::list* RooHistPdf::plotSamplingHint(RooDataHist const& dataHist, // Retrieve position of all bin boundaries const RooAbsBinning* binning = lval->getBinningPtr(nullptr); - std::span boundaries{binning->array(), static_cast(binning->numBoundaries())}; + std::span boundaries{binning->array(), static_cast(binning->numBoundaries())}; - auto hint = new std::list ; - - const double delta = (xhi-xlo)*1e-8 ; - - // Sample points right next to the plot limits - hint->push_back(xlo + delta); - hint->push_back(xhi - delta); - - // Sample points very close to the left and right of the bin boundaries that - // are strictly in between the plot limits. - for (const double x : boundaries) { - if (x - xlo > delta && xhi - x > delta) { - hint->push_back(x - delta); - hint->push_back(x + delta); - } - } - - return hint ; + // Use the helper function from RooCurve to make sure to get sampling hints + // that work with the RooFitPlotting. + return RooCurve::plotSamplingHintForBinBoundaries(boundaries, xlo, xhi); }