From 14e8707f33afb8f8a1ac47f4aace20d8e19bffd6 Mon Sep 17 00:00:00 2001 From: pitkajuh Date: Sat, 6 Jul 2024 12:39:55 +0300 Subject: [PATCH] Fix #884, remove code repetition in ProcessNormalization and AsymPow regarding logKappaForX --- interface/ProcessNormalization.h | 8 ++--- interface/logKappa.h | 23 +++++++++++++++ src/AsymPow.cc | 31 ++++++-------------- src/ProcessNormalization.cc | 50 ++++++++++++-------------------- 4 files changed, 54 insertions(+), 58 deletions(-) create mode 100644 interface/logKappa.h diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index 9a81d6ec498..9e403b53bda 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -33,20 +33,20 @@ class ProcessNormalization : public RooAbsReal { private: // ---- PERSISTENT ---- - double nominalValue_; + double nominalValue_; std::vector logKappa_; // Logarithm of symmetric kappas RooListProxy thetaList_; // List of nuisances for symmetric kappas std::vector > logAsymmKappa_; // Logarithm of asymmetric kappas (low, high) RooListProxy asymmThetaList_; // List of nuisances for asymmetric kappas - RooListProxy otherFactorList_; // Other multiplicative terms + 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 ; + Double_t logKappaForX(double x, const double &kappaHigh, const double &kappaLow) const ; - ClassDefOverride(ProcessNormalization,1) // Process normalization interpolator + ClassDefOverride(ProcessNormalization,1) // Process normalization interpolator }; #endif diff --git a/interface/logKappa.h b/interface/logKappa.h new file mode 100644 index 00000000000..0b66cff1a1f --- /dev/null +++ b/interface/logKappa.h @@ -0,0 +1,23 @@ +#ifndef logKappa_h +#define logKappa_h + +template +T logKappa(const T x, const double &kappaHigh, const double &kappaLow) { + if (fabs(x) >= 0.5) return (x >= 0 ? kappaHigh : -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 + const double logKhi = kappaHigh; + const double logKlo = -kappaLow; + const double avg = 0.5*(logKhi + logKlo), halfdiff = 0.5*(logKhi - logKlo); + const double twox = x+x, twox2 = twox*twox; + const double alpha = 0.125 * twox * (twox2 * (3*twox2 - 10.) + 15.); + return avg + alpha*halfdiff; +} + +#endif diff --git a/src/AsymPow.cc b/src/AsymPow.cc index 04d7a2b0b5d..49106c05b9f 100644 --- a/src/AsymPow.cc +++ b/src/AsymPow.cc @@ -1,4 +1,5 @@ #include "../interface/AsymPow.h" +#include "../interface/logKappa.h" #include #include @@ -6,9 +7,9 @@ AsymPow::AsymPow(const char *name, const char *title, RooAbsReal &kappaLow, RooAbsReal &kappaHigh, RooAbsReal &theta) : RooAbsReal(name,title), - kappaLow_("kappaLow","Base for theta < 0", this, kappaLow), + kappaLow_("kappaLow","Base for theta < 0", this, kappaLow), kappaHigh_("kappaHigh","Base for theta > 0", this, kappaHigh), - theta_("theta", "Exponent (unit gaussian)", this, theta) + theta_("theta", "Exponent (unit gaussian)", this, theta) { } AsymPow::AsymPow(const AsymPow &other, const char *newname) : @@ -20,7 +21,7 @@ AsymPow::AsymPow(const AsymPow &other, const char *newname) : AsymPow::~AsymPow() {} -TObject *AsymPow::clone(const char *newname) const +TObject *AsymPow::clone(const char *newname) const { return new AsymPow(*this,newname); } @@ -31,28 +32,14 @@ Double_t AsymPow::evaluate() const { } Double_t AsymPow::logKappaForX(Double_t x) const { + const double logKhi = log(kappaHigh_); + const double logKlo = -log(kappaLow_); #if 0 // old version with discontinuous derivatives - return (x >= 0 ? log(kappaHigh_) : - log(kappaLow_)); + return (x >= 0 ? logKhi : logKlo); #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; + return logKappa(x, logKhi, logKlo); #endif -} +} ClassImp(AsymPow) diff --git a/src/ProcessNormalization.cc b/src/ProcessNormalization.cc index f61c9b17628..b9ab73e8819 100644 --- a/src/ProcessNormalization.cc +++ b/src/ProcessNormalization.cc @@ -1,4 +1,5 @@ #include "../interface/ProcessNormalization.h" +#include "../interface/logKappa.h" #include #include @@ -7,19 +8,19 @@ ProcessNormalization::ProcessNormalization(const char *name, const char *title, double nominal) : RooAbsReal(name,title), nominalValue_(nominal), - thetaList_("thetaList","List of nuisances for symmetric kappas", this), - asymmThetaList_("asymmThetaList","List of nuisances for asymmetric kappas", this), + thetaList_("thetaList","List of nuisances for symmetric kappas", this), + asymmThetaList_("asymmThetaList","List of nuisances for asymmetric kappas", this), otherFactorList_("otherFactorList","Other multiplicative terms", this) -{ +{ } ProcessNormalization::ProcessNormalization(const char *name, const char *title, RooAbsReal &nominal) : RooAbsReal(name,title), nominalValue_(1.0), - thetaList_("thetaList", "List of nuisances for symmetric kappas", this), - asymmThetaList_("asymmThetaList", "List of nuisances for asymmetric kappas", this), + thetaList_("thetaList", "List of nuisances for symmetric kappas", this), + asymmThetaList_("asymmThetaList", "List of nuisances for asymmetric kappas", this), otherFactorList_("otherFactorList", "Other multiplicative terms", this) -{ +{ otherFactorList_.add(nominal); } @@ -27,9 +28,9 @@ ProcessNormalization::ProcessNormalization(const ProcessNormalization &other, co RooAbsReal(other, newname ? newname : other.GetName()), nominalValue_(other.nominalValue_), logKappa_(other.logKappa_), - thetaList_("thetaList", this, other.thetaList_), + thetaList_("thetaList", this, other.thetaList_), logAsymmKappa_(other.logAsymmKappa_), - asymmThetaList_("asymmThetaList", this, other.asymmThetaList_), + asymmThetaList_("asymmThetaList", this, other.asymmThetaList_), otherFactorList_("otherFactorList", this, other.otherFactorList_) { } @@ -93,7 +94,7 @@ Double_t ProcessNormalization::evaluate() const { const RooAbsReal *theta = asymmThetaListVec_.at(i); const std::pair logKappas = logAsymmKappa_.at(i); double x = theta->getVal(); - logVal += x * logKappaForX(x, logKappas); + logVal += x * logKappaForX(x, logKappas.second, logKappas.first); } } double norm = nominalValue_; @@ -106,43 +107,28 @@ Double_t ProcessNormalization::evaluate() const { 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; -} +Double_t ProcessNormalization::logKappaForX(double x, const double &kappaHigh, const double &kappaLow) const { + return logKappa(x, kappaHigh, kappaLow); +} void ProcessNormalization::dump() const { std::cout << "Dumping ProcessNormalization " << GetName() << " @ " << (void*)this << std::endl; std::cout << "\tnominal value: " << nominalValue_ << std::endl; std::cout << "\tlog-normals (" << logKappa_.size() << "):" << std::endl; for (unsigned int i = 0; i < logKappa_.size(); ++i) { - std::cout << "\t\t kappa = " << exp(logKappa_[i]) << ", logKappa = " << logKappa_[i] << + std::cout << "\t\t kappa = " << exp(logKappa_[i]) << ", logKappa = " << logKappa_[i] << ", theta = " << thetaList_.at(i)->GetName() << " = " << ((RooAbsReal*)thetaList_.at(i))->getVal() << std::endl; } std::cout << "\tasymm log-normals (" << logAsymmKappa_.size() << "):" << std::endl; for (unsigned int i = 0; i < logAsymmKappa_.size(); ++i) { - std::cout << "\t\t kappaLo = " << exp(logAsymmKappa_[i].first) << ", logKappaLo = " << logAsymmKappa_[i].first << - ", kappaHi = " << exp(logAsymmKappa_[i].second) << ", logKappaHi = " << logAsymmKappa_[i].second << + std::cout << "\t\t kappaLo = " << exp(logAsymmKappa_[i].first) << ", logKappaLo = " << logAsymmKappa_[i].first << + ", kappaHi = " << exp(logAsymmKappa_[i].second) << ", logKappaHi = " << logAsymmKappa_[i].second << ", theta = " << asymmThetaList_.at(i)->GetName() << " = " << ((RooAbsReal*)asymmThetaList_.at(i))->getVal() << std::endl; } std::cout << "\tother terms (" << otherFactorList_.getSize() << "):" << std::endl; - for (int i = 0; i < otherFactorList_.getSize(); ++i) { + for (int i = 0; i < otherFactorList_.getSize(); ++i) { std::cout << "\t\t term " << otherFactorList_.at(i)->GetName() << - " (class " << otherFactorList_.at(i)->ClassName() << + " (class " << otherFactorList_.at(i)->ClassName() << "), value = " << ((RooAbsReal*)otherFactorList_.at(i))->getVal() << std::endl; } std::cout << std::endl;