[RF] Add multiplicative 6th order poly + lin. extrapol. (interpCode 6)

Introduce a new interpolation code interpCode=6 for
PiecewiseInterpolation: Multiplicative Polynominal Interpolation with
Linear Extrapolation.

ATLAS Higgs group has a need for interpCode 4 (poly interp with linear
extrap) but in a multiplicative form rather than additive. The new
interpCode 6 adds this option.

I also tried to update the documentation in PiecewiseInterpolation to
explain all the interpcodes (I skipped interpCode 3 since it looks like
its actually identical to interpCode 2 at the moment, and anyway I don't
see it used in the wild at the moment).

roofit/histfactory/src/PiecewiseInterpolation.cxx

@@ -5,17 +5,35 @@
  * more altered or distorted ones, it computes a new shape depending on the value of the nuisance
  * parameters \f$ \alpha_i \f$:
  * \f[
- *   A = \sum_i \mathrm{Interpolate}(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i, \alpha_i).
+ *   A = \mathrm{nominal} + \sum_i I(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i, \alpha_i).
  * \f]
+ * for additive interpolation modes (interpCodes 0, 2, 3, and 4), or:
+ * \f[
+ *   A = \mathrm{nominal}\prod_i I(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i, \alpha_i).
+ * \f]
+ * for multiplicative interpolation modes (interpCodes 1, 5, and 6). The interpCodes determine the function \f$ I \f$ (see table below).
+ * 
+ * Note that a PiecewiseInterpolation with \f$ \mathrm{nominal}=1 \f$, N variations, and a multiplicative interpolation mode is equivalen to N
+ * PiecewiseInterpolations each with a single variation and the same interpolation code, all inside a RooProduct.
+ * 
  * If an \f$ \alpha_i \f$ is zero, the distribution is identical to the nominal distribution, at
  * \f$ \pm 1 \f$ it is identical to the up/down distribution for that specific \f$ i \f$.
  *
+ * PiecewiseInterpolation will behave identically (except for differences in the interpCode assignments) to a FlexibleInterpVar if both its nominal, and high and low variation sets
+ * are all RooRealVar.
+ *
  * The class supports several interpolation methods, which can be selected for each parameter separately
- * using setInterpCode(). The default interpolation code is 4. This performs
- * - \f$ |\alpha | > 1 \f$: Linear extrapolation.
- * - \f$ |\alpha | < 1 \f$: Polynomial interpolation. A sixth-order polynomial is used. Its coefficients
- * are chosen such that function, first, and second derivative at \f$ \alpha \pm 1 \f$ match the values
- * that the extrapolation procedure uses.
+ * using setInterpCode(). The default interpolation code is 0. The table below provides details of the interpCodes:
+
+| **interpCode** | **Name** | **Description** |
+|----------------|----------|-----------------|
+| 0   (default)  | Additive Piecewise Linear | \f$ I = \alpha_i(\mathrm{high}_i - \mathrm{nominal}) \f$ for \f$ \alpha_i>=0 \f$, otherwise \f$ I = \alpha_i(\mathrm{nominal}-\mathrm{low}_i) \f$. Not recommended except if using a symmetric variation, because of discontinuities in derivatives. |
+| 1              | Multiplicative Piecewise Exponential | \f$ I = (\mathrm{high}_i/\mathrm{nominal})^{\alpha_i} \f$ for \f$ \alpha_i>=0 \f$, otherwise \f$ I = (\mathrm{low}_i/\mathrm{nominal})^{-\alpha_i} \f$. |
+| 2              | Additive Quadratic Interp. + Linear Extrap. | Deprecated by interpCode 4. |
+| 4              | Additive Poly Interp. + Linear Extrap. | interpCode 0 for \f$ |\alpha_i|>=1 \f$. 6th-order polynomial for \f$ |\alpha_i|<1 \f$ with matching 0th,1st,2nd derivatives at boundary. Recommended for normalization factors. |
+| 5              | Multiplicative Poly Interp. + Exponential Extrap. | interpCode 1 for \f$ |\alpha_i|>=1 \f$. 6th-order polynomial for \f$ |\alpha_i|<1 \f$ with matching 0th,1st,2nd derivatives at boundary. Recommended for normalization factors. In FlexibleInterpVar this is interpCode=4. |
+| 6              | Multiplicative Poly Interp. + Linear Extrap. | Multiplicative version of interpCode 4. Recommended for normalization factors that must not have roots (i.e. be equal to 0) outside of \f$ |\alpha_i|<1 \f$. |
+
  */
 
 #include "RooStats/HistFactory/PiecewiseInterpolation.h"
@@ -576,7 +594,7 @@ void PiecewiseInterpolation::setAllInterpCodes(int code)
 void PiecewiseInterpolation::setInterpCodeForParam(int iParam, int code)
 {
    RooAbsArg const &param = _paramSet[iParam];
-   if (code < 0 || code > 5) {
+   if (code < 0 || code > 6) {
       coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName()
                             << " with unknown interpolation code " << code << ", keeping current code "
                             << _interpCode[iParam] << std::endl;

roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h

@@ -175,7 +175,7 @@ inline unsigned int getUniformBinning(double low, double high, double val, unsig
    return val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
 }
 
-inline double interpolate1d(double low, double high, double val, unsigned int numBins, double const* vals)
+inline double interpolate1d(double low, double high, double val, unsigned int numBins, double const *vals)
 {
    double binWidth = (high - low) / numBins;
    int idx = val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
@@ -185,9 +185,9 @@ inline double interpolate1d(double low, double high, double val, unsigned int nu
    if (val > low + 0.5 * binWidth && val < high - 0.5 * binWidth) {
       double slope;
       if (val < central) {
-          slope = vals[idx] - vals[idx - 1];
+         slope = vals[idx] - vals[idx - 1];
       } else {
-          slope = vals[idx + 1] - vals[idx];
+         slope = vals[idx + 1] - vals[idx];
       }
       return vals[idx] + slope * (val - central) / binWidth;
    }
@@ -239,10 +239,10 @@ inline double flexibleInterpSingle(unsigned int code, double low, double high, d
       } else {
          return a * std::pow(paramVal, 2) + b * paramVal + c;
       }
-   // According to an old comment in the source code, code 3 was apparently
-   // meant to be a "parabolic version of log-normal", but it never got
-   // implemented. If someone would need it, it could be implemented as doing
-   // code 2 in log space.
+      // According to an old comment in the source code, code 3 was apparently
+      // meant to be a "parabolic version of log-normal", but it never got
+      // implemented. If someone would need it, it could be implemented as doing
+      // code 2 in log space.
    } else if (code == 4) {
       double x = paramVal;
       if (x >= boundary) {
@@ -304,6 +304,24 @@ inline double flexibleInterpSingle(unsigned int code, double low, double high, d
          mod = value;
       }
       return res * (mod - 1.0);
+   } else if (code == 6) {
+      // multiplicative version of code 4 (6th order poly interp with linear extrap)
+
+      double x = paramVal;
+      if (x >= boundary) {
+         return res * (x * (high - nominal) - 1.);
+      } else if (x <= -boundary) {
+         return res * (x * (nominal - low) - 1.);
+      }
+
+      // interpolate 6th degree
+      double t = x / boundary;
+      double eps_plus = high - nominal;
+      double eps_minus = nominal - low;
+      double S = 0.5 * (eps_plus + eps_minus);
+      double A = 0.0625 * (eps_plus - eps_minus);
+
+      return res * (x * (S + t * A * (15 + t * t * (-10 + t * t * 3))) - 1.);
    }
 
    return 0.0;