From 0aa00c940df8bbe6f0d8c035ae6b76305632603f 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/FastTemplate_Old.h | 3 + interface/ProcessNormalization.h | 14 +++-- interface/VerticalInterpHistPdf.h | 3 + src/AsymPow.cc | 30 +--------- src/ProcessNormalization.cc | 94 ++++++++++--------------------- 6 files changed, 48 insertions(+), 99 deletions(-) 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/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;