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

@@ -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	
 };
 

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 ;

interface/ProcessNormalization.h

@@ -28,23 +28,25 @@ 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<double> logKappa_; // Logarithm of symmetric kappas
         RooListProxy thetaList_;        // List of nuisances for symmetric kappas
         std::vector<std::pair<double,double> > logAsymmKappa_; // Logarithm of asymmetric kappas (low, high)
         RooListProxy asymmThetaList_;                           // List of nuisances for asymmetric kappas
         RooListProxy otherFactorList_;     // Other multiplicative terms 
-        std::vector<RooAbsReal *> thetaListVec_; //! Don't serialize me
-        std::vector<RooAbsReal *> asymmThetaListVec_; //! Don't serialize me
-        std::vector<RooAbsReal *> otherFactorListVec_; //! Don't serialize me
-
-        // get the kappa for the appropriate x
-        Double_t logKappaForX(double x, const std::pair<double,double> &logKappas ) const ;
+        mutable std::vector<double> thetaListVec_; //! Don't serialize me
+        mutable std::vector<double> asymmThetaListVec_; //! Don't serialize me
+        mutable std::vector<double> otherFactorListVec_; //! Don't serialize me
+        mutable std::vector<double> logAsymmKappaLow_; //! Don't serialize me
+        mutable std::vector<double> logAsymmKappaHigh_; //! Don't serialize me
 
   ClassDefOverride(ProcessNormalization,1) // Process normalization interpolator 
 };

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<const RooRealVar &>(_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;

src/AsymPow.cc

@@ -1,5 +1,7 @@
 #include "../interface/AsymPow.h"
 
+#include "../interface/CombineMathFuncs.h"
+
 #include <cmath>
 #include <cassert>
 #include <cstdio>
@@ -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)

src/ProcessNormalization.cc

@@ -1,5 +1,7 @@
 #include "../interface/ProcessNormalization.h"
 
+#include "../interface/CombineMathFuncs.h"
+
 #include <cmath>
 #include <cassert>
 #include <cstdio>
@@ -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<RooAbsReal *> & thetaListVec = const_cast<std::vector<RooAbsReal *>&>(thetaListVec_);
-        thetaListVec.reserve(thetaList_.getSize());
-        for (RooAbsArg *a : thetaList_) {
-            thetaListVec.push_back(dynamic_cast<RooAbsReal *>(a));
-        }
-    }
-    if (asymmThetaListVec_.empty()) {
-        std::vector<RooAbsReal *> & asymmThetaListVec = const_cast<std::vector<RooAbsReal *>&>(asymmThetaListVec_);
-        asymmThetaListVec.reserve(asymmThetaList_.getSize());
-        for (RooAbsArg *a : asymmThetaList_) {
-            asymmThetaListVec.push_back(dynamic_cast<RooAbsReal *>(a));
-        }
-    }
-    if (otherFactorListVec_.empty()) {
-        std::vector<RooAbsReal *> & otherFactorListVec = const_cast<std::vector<RooAbsReal *>&>(otherFactorListVec_);
-        otherFactorListVec.reserve(otherFactorList_.getSize());
-        for (RooAbsArg *a : otherFactorList_) {
-            otherFactorListVec.push_back(dynamic_cast<RooAbsReal *>(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<RooAbsReal const&>(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<double,double> 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<RooAbsReal const&>(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<RooAbsReal const&>(otherFactorList_[i]).getVal();
     }
-    return norm;
-}
 
-Double_t ProcessNormalization::logKappaForX(double x, const std::pair<double,double> &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;