From 1ec76088e517f05292c96faab8f60194cde42bd8 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 18 Nov 2024 13:11:55 +0100 Subject: [PATCH] Factor out some mathematical details into free functions This makes it easy to re-use the mathematical functions in generated code. It's part of the `guitargeek/roofit_ad_ichep_2024` branch. --- interface/AsymPow.h | 3 - interface/CombineMathFuncs.h | 152 ++++++++++++++++++++++++++++++ interface/FastTemplate_Old.h | 3 + interface/ProcessNormalization.h | 14 +-- interface/VerticalInterpHistPdf.h | 3 + src/AsymPow.cc | 30 +----- src/ProcessNormalization.cc | 94 ++++++------------ 7 files changed, 200 insertions(+), 99 deletions(-) create mode 100644 interface/CombineMathFuncs.h diff --git a/interface/AsymPow.h b/interface/AsymPow.h index ec38594a2f7..8346bdadba6 100644 --- a/interface/AsymPow.h +++ b/interface/AsymPow.h @@ -37,9 +37,6 @@ class AsymPow : public RooAbsReal { RooRealProxy kappaLow_, kappaHigh_; RooRealProxy theta_; - // get the kappa for the appropriate x - Double_t logKappaForX(Double_t x) const ; - ClassDefOverride(AsymPow,1) // Asymmetric power }; diff --git a/interface/CombineMathFuncs.h b/interface/CombineMathFuncs.h new file mode 100644 index 00000000000..6b461c0d367 --- /dev/null +++ b/interface/CombineMathFuncs.h @@ -0,0 +1,152 @@ +#ifndef CombineMathFuncs_h +#define CombineMathFuncs_h + +#include + +namespace RooFit { +namespace Detail { +namespace MathFuncs { + +inline double smoothStepFunc(double x, double smoothRegion) +{ + if (std::abs(x) >= smoothRegion) + return x > 0 ? +1 : -1; + double xnorm = x / smoothRegion; + double xnorm2 = xnorm * xnorm; + return 0.125 * xnorm * (xnorm2 * (3. * xnorm2 - 10.) + 15); +} + +inline void fastVerticalInterpHistPdf2(int nBins, int nCoefs, double const *coefs, double const *nominal, + double const *binWidth, double const *morphsSum, double const *morphsDiff, + double smoothRegion, double *out) +{ + for (int iBin = 0; iBin < nBins; ++iBin) { + out[iBin] = nominal[iBin]; + } + + double normSum = 0.0; + + for (int iBin = 0; iBin < nBins; ++iBin) { + // apply all morphs one by one + for (int iCoef = 0; iCoef < nCoefs; ++iCoef) { + double const *sum = morphsSum + iCoef * nBins; + double const *diff = morphsDiff + iCoef * nBins; + double x = coefs[iCoef]; + double a = 0.5 * x; + double b = smoothStepFunc(x, smoothRegion); + + out[iBin] += a * (diff[iBin] + b * sum[iBin]); + } + + out[iBin] = std::max(1e-9, out[iBin]); + normSum += out[iBin] * binWidth[iBin]; + } + + if (normSum > 0.0) { + double normSumInv = 1. / normSum; + for (int iBin = 0; iBin < nBins; ++iBin) { + out[iBin] *= normSumInv; + } + } +} + +inline void fastVerticalInterpHistPdf2D2(int nBinsX, int nBinsY, int nCoefs, double const *coefs, + double const *nominal, double const *binWidth, double const *morphsSum, + double const *morphsDiff, double smoothRegion, double *out) +{ + int nBins = nBinsX * nBinsY; + + for (int iBin = 0; iBin < nBins; ++iBin) { + out[iBin] = nominal[iBin]; + } + + for (int iBinX = 0; iBinX < nBinsX; ++iBinX) { + + double normSum = 0.0; + + for (int iBinY = 0; iBinY < nBinsY; ++iBinY) { + int iBin = iBinY + nBinsY * iBinX; + // apply all morphs one by one + for (int iCoef = 0; iCoef < nCoefs; ++iCoef) { + double const *sum = morphsSum + iCoef * nBinsY; + double const *diff = morphsDiff + iCoef * nBinsY; + double x = coefs[iCoef]; + double a = 0.5 * x; + double b = smoothStepFunc(x, smoothRegion); + + out[iBin] += a * (diff[iBin] + b * sum[iBin]); + } + + out[iBin] = std::max(1e-9, out[iBin]); + normSum += out[iBin] * binWidth[iBin]; + } + + if (normSum > 0.0) { + double normSumInv = 1. / normSum; + for (int iBinY = 0; iBinY < nBinsY; ++iBinY) { + int iBin = iBinY + nBinsY * iBinX; + out[iBin] *= normSumInv; + } + } + } +} + +inline double logKappaForX(double theta, double logKappaLow, double logKappaHigh) +{ + double logKappa = 0.0; + + if (std::abs(theta) >= 0.5) { + logKappa = theta >= 0 ? logKappaHigh : -logKappaLow; + } else { + // interpolate between log(kappaHigh) and -log(kappaLow) + // logKappa(x) = avg + halfdiff * h(2x) + // where h(x) is the 3th order polynomial + // h(x) = (3 x^5 - 10 x^3 + 15 x)/8; + // chosen so that h(x) satisfies the following: + // h (+/-1) = +/-1 + // h'(+/-1) = 0 + // h"(+/-1) = 0 + double logKhi = logKappaHigh; + double logKlo = -logKappaLow; + double avg = 0.5 * (logKhi + logKlo); + double halfdiff = 0.5 * (logKhi - logKlo); + double twox = theta + theta; + double twox2 = twox * twox; + double alpha = 0.125 * twox * (twox2 * (3 * twox2 - 10.) + 15.); + logKappa = avg + alpha * halfdiff; + } + + return logKappa; +} + +inline double asymPow(double theta, double kappaLow, double kappaHigh) +{ + return std::exp(logKappaForX(theta, std::log(kappaLow), std::log(kappaHigh)) * theta); +} + +inline double processNormalization(double nominalValue, std::size_t nThetas, std::size_t nAsymmThetas, + std::size_t nOtherFactors, double const *thetas, double const *logKappas, + double const *asymmThetas, double const *asymmLogKappasLow, + double const *asymmLogKappasHigh, double const *otherFactors) +{ + double logVal = 0.0; + for (std::size_t i = 0; i < nThetas; i++) { + logVal += thetas[i] * logKappas[i]; + } + for (std::size_t i = 0; i < nAsymmThetas; i++) { + double x = asymmThetas[i]; + logVal += x * logKappaForX(x, asymmLogKappasLow[i], asymmLogKappasHigh[i]); + } + double norm = nominalValue; + norm *= std::exp(logVal); + for (std::size_t i = 0; i < nOtherFactors; i++) { + norm *= otherFactors[i]; + } + return norm; +} + +} // namespace MathFuncs +} // namespace Detail +} // namespace RooFit + +#endif diff --git a/interface/FastTemplate_Old.h b/interface/FastTemplate_Old.h index a272a559f58..5461d8f2109 100644 --- a/interface/FastTemplate_Old.h +++ b/interface/FastTemplate_Old.h @@ -169,6 +169,9 @@ class FastHisto2D : public FastTemplate { void Dump() const ; + const T & GetBinContent(int bin) const { return values_[bin]; } + const T & GetWidth(unsigned int i) const { return binWidths_[i]; } + T GetMaxOnXY() const ; T GetMaxOnX(const T &y) const ; T GetMaxOnY(const T &x) const ; diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index 9a81d6ec498..24916cc3792 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -28,10 +28,13 @@ class ProcessNormalization : public RooAbsReal { void addAsymmLogNormal(double kappaLo, double kappaHi, RooAbsReal &theta) ; void addOtherFactor(RooAbsReal &factor) ; void dump() const ; + protected: Double_t evaluate() const override; private: + void fillAsymmKappaVecs() const; + // ---- PERSISTENT ---- double nominalValue_; std::vector logKappa_; // Logarithm of symmetric kappas @@ -39,12 +42,11 @@ class ProcessNormalization : public RooAbsReal { std::vector > logAsymmKappa_; // Logarithm of asymmetric kappas (low, high) RooListProxy asymmThetaList_; // List of nuisances for asymmetric kappas RooListProxy otherFactorList_; // Other multiplicative terms - std::vector thetaListVec_; //! Don't serialize me - std::vector asymmThetaListVec_; //! Don't serialize me - std::vector otherFactorListVec_; //! Don't serialize me - - // get the kappa for the appropriate x - Double_t logKappaForX(double x, const std::pair &logKappas ) const ; + mutable std::vector thetaListVec_; //! Don't serialize me + mutable std::vector asymmThetaListVec_; //! Don't serialize me + mutable std::vector otherFactorListVec_; //! Don't serialize me + mutable std::vector logAsymmKappaLow_; //! Don't serialize me + mutable std::vector logAsymmKappaHigh_; //! Don't serialize me ClassDefOverride(ProcessNormalization,1) // Process normalization interpolator }; diff --git a/interface/VerticalInterpHistPdf.h b/interface/VerticalInterpHistPdf.h index c07a3aa5e2e..f1172b6f6d0 100644 --- a/interface/VerticalInterpHistPdf.h +++ b/interface/VerticalInterpHistPdf.h @@ -304,6 +304,8 @@ class FastVerticalInterpHistPdf2 : public FastVerticalInterpHistPdf2Base { TObject* clone(const char* newname) const override { return new FastVerticalInterpHistPdf2(*this,newname) ; } ~FastVerticalInterpHistPdf2() override {} + const RooRealVar & x() const { return dynamic_cast(_x.arg()); } + void setActiveBins(unsigned int bins) override ; Double_t evaluate() const override ; @@ -366,6 +368,7 @@ class FastVerticalInterpHistPdf2D2 : public FastVerticalInterpHistPdf2Base { Double_t maxVal(Int_t code) const override ; Double_t evaluate() const override ; + protected: RooRealProxy _x, _y; bool _conditional; diff --git a/src/AsymPow.cc b/src/AsymPow.cc index 04d7a2b0b5d..fb7ba951c13 100644 --- a/src/AsymPow.cc +++ b/src/AsymPow.cc @@ -1,5 +1,7 @@ #include "../interface/AsymPow.h" +#include "../interface/CombineMathFuncs.h" + #include #include #include @@ -26,33 +28,7 @@ TObject *AsymPow::clone(const char *newname) const } Double_t AsymPow::evaluate() const { - Double_t x = theta_; - return exp(logKappaForX(x) * x); + return RooFit::Detail::MathFuncs::asymPow(theta_, kappaLow_, kappaHigh_); } -Double_t AsymPow::logKappaForX(Double_t x) const { -#if 0 - // old version with discontinuous derivatives - return (x >= 0 ? log(kappaHigh_) : - log(kappaLow_)); -#else - if (fabs(x) >= 0.5) return (x >= 0 ? log(kappaHigh_) : - log(kappaLow_)); - // interpolate between log(kappaHigh) and -log(kappaLow) - // logKappa(x) = avg + halfdiff * h(2x) - // where h(x) is the 3th order polynomial - // h(x) = (3 x^5 - 10 x^3 + 15 x)/8; - // chosen so that h(x) satisfies the following: - // h (+/-1) = +/-1 - // h'(+/-1) = 0 - // h"(+/-1) = 0 - double logKhi = log(kappaHigh_); - double logKlo = -log(kappaLow_); - double avg = 0.5*(logKhi + logKlo), halfdiff = 0.5*(logKhi - logKlo); - double twox = x+x, twox2 = twox*twox; - double alpha = 0.125 * twox * (twox2 * (3*twox2 - 10.) + 15.); - double ret = avg + alpha*halfdiff; - //assert(alpha >= -1 && alpha <= 1 && "Something is wrong in the interpolation"); - return ret; -#endif -} - ClassImp(AsymPow) diff --git a/src/ProcessNormalization.cc b/src/ProcessNormalization.cc index f61c9b17628..5baa229b3a1 100644 --- a/src/ProcessNormalization.cc +++ b/src/ProcessNormalization.cc @@ -1,5 +1,7 @@ #include "../interface/ProcessNormalization.h" +#include "../interface/CombineMathFuncs.h" + #include #include #include @@ -56,74 +58,40 @@ void ProcessNormalization::addOtherFactor(RooAbsReal &factor) { otherFactorList_.add(factor); } -Double_t ProcessNormalization::evaluate() const { - double logVal = 0.0; - if (thetaListVec_.empty()) { - std::vector & thetaListVec = const_cast&>(thetaListVec_); - thetaListVec.reserve(thetaList_.getSize()); - for (RooAbsArg *a : thetaList_) { - thetaListVec.push_back(dynamic_cast(a)); - } - } - if (asymmThetaListVec_.empty()) { - std::vector & asymmThetaListVec = const_cast&>(asymmThetaListVec_); - asymmThetaListVec.reserve(asymmThetaList_.getSize()); - for (RooAbsArg *a : asymmThetaList_) { - asymmThetaListVec.push_back(dynamic_cast(a)); - } - } - if (otherFactorListVec_.empty()) { - std::vector & otherFactorListVec = const_cast&>(otherFactorListVec_); - otherFactorListVec.reserve(otherFactorList_.getSize()); - for (RooAbsArg *a : otherFactorList_) { - otherFactorListVec.push_back(dynamic_cast(a)); - } +void ProcessNormalization::fillAsymmKappaVecs() const +{ + if (logAsymmKappaLow_.size() != logAsymmKappa_.size()) { + logAsymmKappaLow_.reserve(logAsymmKappa_.size()); + logAsymmKappaLow_.reserve(logAsymmKappa_.size()); + for (auto [lo, hi] : logAsymmKappa_) { + logAsymmKappaLow_.push_back(lo); + logAsymmKappaHigh_.push_back(hi); + } } - if (!logKappa_.empty()) { - assert(logKappa_.size()==thetaListVec_.size()); - for(unsigned int i=0; i < thetaListVec_.size() ; i++){ - const RooAbsReal *theta = thetaListVec_.at(i); - const double logKappa = logKappa_.at(i); - logVal += theta->getVal() * (logKappa); - } +} + +Double_t ProcessNormalization::evaluate() const +{ + thetaListVec_.resize(thetaList_.size()); + asymmThetaListVec_.resize(asymmThetaList_.size()); + otherFactorListVec_.resize(otherFactorList_.size()); + for (std::size_t i = 0; i < thetaList_.size(); ++i) { + thetaListVec_[i] = static_cast(thetaList_[i]).getVal(); } - if (!logAsymmKappa_.empty()) { - assert(logAsymmKappa_.size()==asymmThetaListVec_.size()); - for( unsigned int i=0; i < asymmThetaListVec_.size(); i++){ - const RooAbsReal *theta = asymmThetaListVec_.at(i); - const std::pair logKappas = logAsymmKappa_.at(i); - double x = theta->getVal(); - logVal += x * logKappaForX(x, logKappas); - } + for (std::size_t i = 0; i < asymmThetaList_.size(); ++i) { + asymmThetaListVec_[i] = static_cast(asymmThetaList_[i]).getVal(); } - double norm = nominalValue_; - if (logVal) norm *= std::exp(logVal); - if (otherFactorList_.getSize()) { - for (const RooAbsReal *fact :otherFactorListVec_){ - norm *= fact->getVal(); - } + for (std::size_t i = 0; i < otherFactorList_.size(); ++i) { + otherFactorListVec_[i] = static_cast(otherFactorList_[i]).getVal(); } - return norm; -} -Double_t ProcessNormalization::logKappaForX(double x, const std::pair &logKappas) const { - if (fabs(x) >= 0.5) return (x >= 0 ? logKappas.second : -logKappas.first); - // interpolate between log(kappaHigh) and -log(kappaLow) - // logKappa(x) = avg + halfdiff * h(2x) - // where h(x) is the 3th order polynomial - // h(x) = (3 x^5 - 10 x^3 + 15 x)/8; - // chosen so that h(x) satisfies the following: - // h (+/-1) = +/-1 - // h'(+/-1) = 0 - // h"(+/-1) = 0 - double logKhi = logKappas.second; - double logKlo = -logKappas.first; - double avg = 0.5*(logKhi + logKlo), halfdiff = 0.5*(logKhi - logKlo); - double twox = x+x, twox2 = twox*twox; - double alpha = 0.125 * twox * (twox2 * (3*twox2 - 10.) + 15.); - double ret = avg + alpha*halfdiff; - return ret; -} + fillAsymmKappaVecs(); + return RooFit::Detail::MathFuncs::processNormalization(nominalValue_, + thetaList_.size(), asymmThetaList_.size(), otherFactorList_.size(), + thetaListVec_.data(), logKappa_.data(), asymmThetaListVec_.data(), + logAsymmKappaLow_.data(), logAsymmKappaHigh_.data(), + otherFactorListVec_.data()); +} void ProcessNormalization::dump() const { std::cout << "Dumping ProcessNormalization " << GetName() << " @ " << (void*)this << std::endl;