Fix #884, remove code repetition in ProcessNormalization and AsymPow

interface/ProcessNormalization.h

@@ -33,20 +33,20 @@ class ProcessNormalization : public RooAbsReal {
 
     private:
         // ---- PERSISTENT ----
-        double nominalValue_;                         
+        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 
+        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 ;
+        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

interface/logKappa.h

@@ -0,0 +1,23 @@
+#ifndef logKappa_h
+#define logKappa_h
+
+template<typename T>
+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

src/AsymPow.cc

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

src/ProcessNormalization.cc

@@ -1,4 +1,5 @@
 #include "../interface/ProcessNormalization.h"
+#include "../interface/logKappa.h"
 
 #include <cmath>
 #include <cassert>
@@ -7,29 +8,29 @@
 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);
 }
 
 ProcessNormalization::ProcessNormalization(const ProcessNormalization &other, const char *newname) :
         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<double,double> 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<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;
-} 
+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;