InterpolationTables/interface/interpolate.h

0001 #ifndef NPSTAT_INTERPOLATE_HH_
0002 #define NPSTAT_INTERPOLATE_HH_
0003
0004 /*!
0005 // \file interpolate.h
0006 //
0007 // \brief Low-order polynomial interpolation (up to cubic)
0008 //        for equidistant coordinates
0009 //
0010 // Author: I. Volobouev
0011 //
0012 // October 2009
0013 */
0014
0015 #include "JetMETCorrections/InterpolationTables/interface/ProperDblFromCmpl.h"
0016
0017 namespace npstat {
0018   /**
0019     // Linear interpolation. Assumes that
0020     // the function values are given at 0 and 1.
0021     */
0022   template <typename T>
0023   inline T interpolate_linear(const double x, const T& f0, const T& f1) {
0024     const typename ProperDblFromCmpl<T>::type dx = 1.0 - x;
0025     return f0 * dx + f1 * static_cast<typename ProperDblFromCmpl<T>::type>(x);
0026   }
0027
0028   /**
0029     // Quadratic interpolation. Assume that
0030     // the function values are given at 0, 1, and 2.
0031     */
0032   template <typename T>
0033   inline T interpolate_quadratic(const double x, const T& f0, const T& f1, const T& f2) {
0034     static const typename ProperDblFromCmpl<T>::type two = 2.0;
0035     const typename ProperDblFromCmpl<T>::type dx = x - 1.0;
0036     return f1 + ((f2 - f0) / two + ((f2 - f1) + (f0 - f1)) * (dx / two)) * dx;
0037   }
0038
0039   /**
0040     // Cubic interpolation. Assume that
0041     // the function values are given at 0, 1, 2, and 3.
0042     */
0043   template <typename T>
0044   inline T interpolate_cubic(const double x, const T& f0, const T& f1, const T& f2, const T& f3) {
0045     return interpolate_linear(
0046         x * (3.0 - x) / 2.0, interpolate_linear(x / 3.0, f0, f3), interpolate_linear(x - 1.0, f1, f2));
0047   }
0048
0049   //@{
0050   /**
0051     // Get the coefficients of the interpolating polynomial.
0052     // The interpolated function values are provided at 0, 1, ...
0053     // The return value of the function is the number of
0054     // coefficients (i.e., the polynomial degree plus one).
0055     // On exit, the coefficients are placed into the "buffer"
0056     // array in the order of increasing monomial degree.
0057     // The length of the provided buffer must be sufficient
0058     // to hold all these coefficients.
0059     */
0060   template <typename T>
0061   unsigned interpolation_coefficients(T* buffer, unsigned bufLen, const T& f0, const T& f1);
0062   template <typename T>
0063   unsigned interpolation_coefficients(T* buffer, unsigned bufLen, const T& f0, const T& f1, const T& f2);
0064   template <typename T>
0065   unsigned interpolation_coefficients(T* buffer, unsigned bufLen, const T& f0, const T& f1, const T& f2, const T& f3);
0066   //@}
0067 }  // namespace npstat
0068
0069 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0070
0071 namespace npstat {
0072   template <typename T>
0073   unsigned interpolation_coefficients(T* buffer, const unsigned bufLen, const T& f0, const T& f1) {
0074     if (bufLen <= 1U)
0075       throw npstat::NpstatInvalidArgument(
0076           "In npstat::interpolation_coefficients: "
0077           "insufficient length of the output buffer");
0078     buffer[0] = f0;
0079     buffer[1] = f1 - f0;
0080     return 2U;
0081   }
0082
0083   template <typename T>
0084   unsigned interpolation_coefficients(T* buffer, const unsigned bufLen, const T& f0, const T& f1, const T& f2) {
0085     if (bufLen <= 2U)
0086       throw npstat::NpstatInvalidArgument(
0087           "In npstat::interpolation_coefficients: "
0088           "insufficient length of the output buffer");
0089     buffer[0] = f0;
0090     buffer[1] = static_cast<T>((f1 - f2 + 3 * (f1 - f0)) / 2.0);
0091     buffer[2] = static_cast<T>(((f0 - f1) + (f2 - f1)) / 2.0);
0092     return 3U;
0093   }
0094
0095   template <typename T>
0096   unsigned interpolation_coefficients(
0097       T* buffer, const unsigned bufLen, const T& f0, const T& f1, const T& f2, const T& f3) {
0098     if (bufLen <= 3U)
0099       throw npstat::NpstatInvalidArgument(
0100           "In npstat::interpolation_coefficients: "
0101           "insufficient length of the output buffer");
0102     buffer[0] = f0;
0103     buffer[1] = static_cast<T>((11 * (f1 - f0) + 7 * (f1 - f2) + 2 * (f3 - f2)) / 6.0);
0104     buffer[2] = static_cast<T>((2 * (f0 - f1) + 3 * (f2 - f1) + (f2 - f3)) / 2.0);
0105     buffer[3] = static_cast<T>(((f3 - f0) + 3 * (f1 - f2)) / 6.0);
0106     return 4U;
0107   }
0108 }  // namespace npstat
0109
0110 #endif  // NPSTAT_INTERPOLATE_HH_