Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:19:21

0001 #ifndef NPSTAT_INTERPOLATEHISTOND_HH_
0002 #define NPSTAT_INTERPOLATEHISTOND_HH_
0003 
0004 /*!
0005 // \file interpolateHistoND.h
0006 //
0007 // \brief Interpolate histogram contents
0008 //
0009 // Functions which interpolate histogram contents are not included
0010 // into the HistoND template itself because we do not always want to
0011 // create histograms using bin types which can be multiplied by doubles
0012 // (also, results of such a multiplication have to be automatically
0013 // converted back to the same type).
0014 //
0015 // The implementations work by invoking "interpolate1" or "interpolate3"
0016 // ArrayND methods on the histogram bin contents after an appropriate
0017 // coordinate transformation.
0018 //
0019 // Author: I. Volobouev
0020 //
0021 // November 2011
0022 */
0023 
0024 #include "JetMETCorrections/InterpolationTables/interface/HistoND.h"
0025 
0026 namespace npstat {
0027   /**
0028     // The interpolation degree in this method can be set to 0, 1, or 3
0029     // which results, respectively, in closest bin lookup, multilinear
0030     // interpolation, or multicubic interpolation. Value of the closest
0031     // bin inside the histogram range is used if some coordinate is outside
0032     // of the corresponding axis limits.
0033     */
0034   template <typename Float, class Axis>
0035   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0036                            const double* coords,
0037                            unsigned coordsDim,
0038                            unsigned interpolationDegree);
0039   //@{
0040   /**
0041     // Convenience function for interpolating histograms, with
0042     // an explicit coordinate argument for each histogram dimension
0043     */
0044   template <typename Float, class Axis>
0045   Float interpolateHistoND(const HistoND<Float, Axis>& histo, double x0, unsigned interpolationDegree);
0046 
0047   template <typename Float, class Axis>
0048   Float interpolateHistoND(const HistoND<Float, Axis>& histo, double x0, double x1, unsigned interpolationDegree);
0049 
0050   template <typename Float, class Axis>
0051   Float interpolateHistoND(
0052       const HistoND<Float, Axis>& histo, double x0, double x1, double x2, unsigned interpolationDegree);
0053 
0054   template <typename Float, class Axis>
0055   Float interpolateHistoND(
0056       const HistoND<Float, Axis>& histo, double x0, double x1, double x2, double x3, unsigned interpolationDegree);
0057 
0058   template <typename Float, class Axis>
0059   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0060                            double x0,
0061                            double x1,
0062                            double x2,
0063                            double x3,
0064                            double x4,
0065                            unsigned interpolationDegree);
0066 
0067   template <typename Float, class Axis>
0068   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0069                            double x0,
0070                            double x1,
0071                            double x2,
0072                            double x3,
0073                            double x4,
0074                            double x5,
0075                            unsigned interpolationDegree);
0076 
0077   template <typename Float, class Axis>
0078   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0079                            double x0,
0080                            double x1,
0081                            double x2,
0082                            double x3,
0083                            double x4,
0084                            double x5,
0085                            double x6,
0086                            unsigned interpolationDegree);
0087 
0088   template <typename Float, class Axis>
0089   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0090                            double x0,
0091                            double x1,
0092                            double x2,
0093                            double x3,
0094                            double x4,
0095                            double x5,
0096                            double x6,
0097                            double x7,
0098                            unsigned interpolationDegree);
0099 
0100   template <typename Float, class Axis>
0101   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0102                            double x0,
0103                            double x1,
0104                            double x2,
0105                            double x3,
0106                            double x4,
0107                            double x5,
0108                            double x6,
0109                            double x7,
0110                            double x8,
0111                            unsigned interpolationDegree);
0112 
0113   template <typename Float, class Axis>
0114   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0115                            double x0,
0116                            double x1,
0117                            double x2,
0118                            double x3,
0119                            double x4,
0120                            double x5,
0121                            double x6,
0122                            double x7,
0123                            double x8,
0124                            double x9,
0125                            unsigned interpolationDegree);
0126   //@}
0127 }  // namespace npstat
0128 
0129 #include <cassert>
0130 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0131 
0132 namespace npstat {
0133   namespace Private {
0134     template <typename Float, class Axis>
0135     void iHND_checkArgs(const HistoND<Float, Axis>& histo, const unsigned xDim, const unsigned interpolationDegree) {
0136       if (xDim != histo.dim())
0137         throw npstat::NpstatInvalidArgument(
0138             "In npstat::interpolateHistoND: incompatible "
0139             "dimensionality of input coordinates");
0140       if (xDim == 0U)
0141         throw npstat::NpstatInvalidArgument(
0142             "In npstat::interpolateHistoND: can not interpolate "
0143             "zero-dimensional histograms");
0144       if (!(interpolationDegree == 0U || interpolationDegree == 1U || interpolationDegree == 3U))
0145         throw npstat::NpstatInvalidArgument(
0146             "In npstat::interpolateHistoND: "
0147             "unsupported interpolation degree");
0148       if (interpolationDegree == 3U && !histo.isUniformlyBinned())
0149         throw npstat::NpstatInvalidArgument(
0150             "In npstat::interpolateHistoND: unsupported "
0151             "interpolation degree for non-uniform binning");
0152     }
0153   }  // namespace Private
0154 
0155   template <typename Float, class Axis>
0156   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0157                            const double* x,
0158                            const unsigned xDim,
0159                            const unsigned interpolationDegree) {
0160     Private::iHND_checkArgs(histo, xDim, interpolationDegree);
0161     assert(x);
0162     const Axis* ax = &histo.axes()[0];
0163     double coords[CHAR_BIT * sizeof(unsigned long)];
0164     for (unsigned i = 0; i < xDim; ++i)
0165       coords[i] = ax[i].fltBinNumber(x[i], false);
0166     const ArrayND<Float>& bins(histo.binContents());
0167     switch (interpolationDegree) {
0168       case 1U:
0169         return bins.interpolate1(coords, xDim);
0170       case 3U:
0171         return bins.interpolate3(coords, xDim);
0172       default:
0173         return bins.closest(coords, xDim);
0174     }
0175   }
0176 
0177   template <typename Float, class Axis>
0178   Float interpolateHistoND(const HistoND<Float, Axis>& histo, const double x0, const unsigned interpolationDegree) {
0179     const unsigned expDim = 1U;
0180     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0181     const double coords = histo.axis(0).fltBinNumber(x0, false);
0182     const ArrayND<Float>& bins(histo.binContents());
0183     switch (interpolationDegree) {
0184       case 1U:
0185         return bins.interpolate1(&coords, expDim);
0186       case 3U:
0187         return bins.interpolate3(&coords, expDim);
0188       default:
0189         return bins.closest(&coords, expDim);
0190     }
0191   }
0192 
0193   template <typename Float, class Axis>
0194   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0195                            const double x0,
0196                            const double x1,
0197                            const unsigned interpolationDegree) {
0198     const unsigned expDim = 2U;
0199     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0200     const Axis* ax = &histo.axes()[0];
0201     double coords[expDim];
0202     coords[0] = ax[0].fltBinNumber(x0, false);
0203     coords[1] = ax[1].fltBinNumber(x1, false);
0204     const ArrayND<Float>& bins(histo.binContents());
0205     switch (interpolationDegree) {
0206       case 1U:
0207         return bins.interpolate1(coords, expDim);
0208       case 3U:
0209         return bins.interpolate3(coords, expDim);
0210       default:
0211         return bins.closest(coords, expDim);
0212     }
0213   }
0214 
0215   template <typename Float, class Axis>
0216   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0217                            const double x0,
0218                            const double x1,
0219                            const double x2,
0220                            const unsigned interpolationDegree) {
0221     const unsigned expDim = 3U;
0222     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0223     const Axis* ax = &histo.axes()[0];
0224     double coords[expDim];
0225     coords[0] = ax[0].fltBinNumber(x0, false);
0226     coords[1] = ax[1].fltBinNumber(x1, false);
0227     coords[2] = ax[2].fltBinNumber(x2, false);
0228     const ArrayND<Float>& bins(histo.binContents());
0229     switch (interpolationDegree) {
0230       case 1U:
0231         return bins.interpolate1(coords, expDim);
0232       case 3U:
0233         return bins.interpolate3(coords, expDim);
0234       default:
0235         return bins.closest(coords, expDim);
0236     }
0237   }
0238 
0239   template <typename Float, class Axis>
0240   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0241                            const double x0,
0242                            const double x1,
0243                            const double x2,
0244                            const double x3,
0245                            const unsigned interpolationDegree) {
0246     const unsigned expDim = 4U;
0247     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0248     const Axis* ax = &histo.axes()[0];
0249     double coords[expDim];
0250     coords[0] = ax[0].fltBinNumber(x0, false);
0251     coords[1] = ax[1].fltBinNumber(x1, false);
0252     coords[2] = ax[2].fltBinNumber(x2, false);
0253     coords[3] = ax[3].fltBinNumber(x3, false);
0254     const ArrayND<Float>& bins(histo.binContents());
0255     switch (interpolationDegree) {
0256       case 1U:
0257         return bins.interpolate1(coords, expDim);
0258       case 3U:
0259         return bins.interpolate3(coords, expDim);
0260       default:
0261         return bins.closest(coords, expDim);
0262     }
0263   }
0264 
0265   template <typename Float, class Axis>
0266   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0267                            const double x0,
0268                            const double x1,
0269                            const double x2,
0270                            const double x3,
0271                            const double x4,
0272                            const unsigned interpolationDegree) {
0273     const unsigned expDim = 5U;
0274     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0275     const Axis* ax = &histo.axes()[0];
0276     double coords[expDim];
0277     coords[0] = ax[0].fltBinNumber(x0, false);
0278     coords[1] = ax[1].fltBinNumber(x1, false);
0279     coords[2] = ax[2].fltBinNumber(x2, false);
0280     coords[3] = ax[3].fltBinNumber(x3, false);
0281     coords[4] = ax[4].fltBinNumber(x4, false);
0282     const ArrayND<Float>& bins(histo.binContents());
0283     switch (interpolationDegree) {
0284       case 1U:
0285         return bins.interpolate1(coords, expDim);
0286       case 3U:
0287         return bins.interpolate3(coords, expDim);
0288       default:
0289         return bins.closest(coords, expDim);
0290     }
0291   }
0292 
0293   template <typename Float, class Axis>
0294   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0295                            const double x0,
0296                            const double x1,
0297                            const double x2,
0298                            const double x3,
0299                            const double x4,
0300                            const double x5,
0301                            const unsigned interpolationDegree) {
0302     const unsigned expDim = 6U;
0303     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0304     const Axis* ax = &histo.axes()[0];
0305     double coords[expDim];
0306     coords[0] = ax[0].fltBinNumber(x0, false);
0307     coords[1] = ax[1].fltBinNumber(x1, false);
0308     coords[2] = ax[2].fltBinNumber(x2, false);
0309     coords[3] = ax[3].fltBinNumber(x3, false);
0310     coords[4] = ax[4].fltBinNumber(x4, false);
0311     coords[5] = ax[5].fltBinNumber(x5, false);
0312     const ArrayND<Float>& bins(histo.binContents());
0313     switch (interpolationDegree) {
0314       case 1U:
0315         return bins.interpolate1(coords, expDim);
0316       case 3U:
0317         return bins.interpolate3(coords, expDim);
0318       default:
0319         return bins.closest(coords, expDim);
0320     }
0321   }
0322 
0323   template <typename Float, class Axis>
0324   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0325                            const double x0,
0326                            const double x1,
0327                            const double x2,
0328                            const double x3,
0329                            const double x4,
0330                            const double x5,
0331                            const double x6,
0332                            const unsigned interpolationDegree) {
0333     const unsigned expDim = 7U;
0334     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0335     const Axis* ax = &histo.axes()[0];
0336     double coords[expDim];
0337     coords[0] = ax[0].fltBinNumber(x0, false);
0338     coords[1] = ax[1].fltBinNumber(x1, false);
0339     coords[2] = ax[2].fltBinNumber(x2, false);
0340     coords[3] = ax[3].fltBinNumber(x3, false);
0341     coords[4] = ax[4].fltBinNumber(x4, false);
0342     coords[5] = ax[5].fltBinNumber(x5, false);
0343     coords[6] = ax[6].fltBinNumber(x6, false);
0344     const ArrayND<Float>& bins(histo.binContents());
0345     switch (interpolationDegree) {
0346       case 1U:
0347         return bins.interpolate1(coords, expDim);
0348       case 3U:
0349         return bins.interpolate3(coords, expDim);
0350       default:
0351         return bins.closest(coords, expDim);
0352     }
0353   }
0354 
0355   template <typename Float, class Axis>
0356   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0357                            const double x0,
0358                            const double x1,
0359                            const double x2,
0360                            const double x3,
0361                            const double x4,
0362                            const double x5,
0363                            const double x6,
0364                            const double x7,
0365                            const unsigned interpolationDegree) {
0366     const unsigned expDim = 8U;
0367     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0368     const Axis* ax = &histo.axes()[0];
0369     double coords[expDim];
0370     coords[0] = ax[0].fltBinNumber(x0, false);
0371     coords[1] = ax[1].fltBinNumber(x1, false);
0372     coords[2] = ax[2].fltBinNumber(x2, false);
0373     coords[3] = ax[3].fltBinNumber(x3, false);
0374     coords[4] = ax[4].fltBinNumber(x4, false);
0375     coords[5] = ax[5].fltBinNumber(x5, false);
0376     coords[6] = ax[6].fltBinNumber(x6, false);
0377     coords[7] = ax[7].fltBinNumber(x7, false);
0378     const ArrayND<Float>& bins(histo.binContents());
0379     switch (interpolationDegree) {
0380       case 1U:
0381         return bins.interpolate1(coords, expDim);
0382       case 3U:
0383         return bins.interpolate3(coords, expDim);
0384       default:
0385         return bins.closest(coords, expDim);
0386     }
0387   }
0388 
0389   template <typename Float, class Axis>
0390   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0391                            const double x0,
0392                            const double x1,
0393                            const double x2,
0394                            const double x3,
0395                            const double x4,
0396                            const double x5,
0397                            const double x6,
0398                            const double x7,
0399                            const double x8,
0400                            const unsigned interpolationDegree) {
0401     const unsigned expDim = 9U;
0402     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0403     const Axis* ax = &histo.axes()[0];
0404     double coords[expDim];
0405     coords[0] = ax[0].fltBinNumber(x0, false);
0406     coords[1] = ax[1].fltBinNumber(x1, false);
0407     coords[2] = ax[2].fltBinNumber(x2, false);
0408     coords[3] = ax[3].fltBinNumber(x3, false);
0409     coords[4] = ax[4].fltBinNumber(x4, false);
0410     coords[5] = ax[5].fltBinNumber(x5, false);
0411     coords[6] = ax[6].fltBinNumber(x6, false);
0412     coords[7] = ax[7].fltBinNumber(x7, false);
0413     coords[8] = ax[8].fltBinNumber(x8, false);
0414     const ArrayND<Float>& bins(histo.binContents());
0415     switch (interpolationDegree) {
0416       case 1U:
0417         return bins.interpolate1(coords, expDim);
0418       case 3U:
0419         return bins.interpolate3(coords, expDim);
0420       default:
0421         return bins.closest(coords, expDim);
0422     }
0423   }
0424 
0425   template <typename Float, class Axis>
0426   Float interpolateHistoND(const HistoND<Float, Axis>& histo,
0427                            const double x0,
0428                            const double x1,
0429                            const double x2,
0430                            const double x3,
0431                            const double x4,
0432                            const double x5,
0433                            const double x6,
0434                            const double x7,
0435                            const double x8,
0436                            const double x9,
0437                            const unsigned interpolationDegree) {
0438     const unsigned expDim = 10U;
0439     Private::iHND_checkArgs(histo, expDim, interpolationDegree);
0440     const Axis* ax = &histo.axes()[0];
0441     double coords[expDim];
0442     coords[0] = ax[0].fltBinNumber(x0, false);
0443     coords[1] = ax[1].fltBinNumber(x1, false);
0444     coords[2] = ax[2].fltBinNumber(x2, false);
0445     coords[3] = ax[3].fltBinNumber(x3, false);
0446     coords[4] = ax[4].fltBinNumber(x4, false);
0447     coords[5] = ax[5].fltBinNumber(x5, false);
0448     coords[6] = ax[6].fltBinNumber(x6, false);
0449     coords[7] = ax[7].fltBinNumber(x7, false);
0450     coords[8] = ax[8].fltBinNumber(x8, false);
0451     coords[9] = ax[9].fltBinNumber(x9, false);
0452     const ArrayND<Float>& bins(histo.binContents());
0453     switch (interpolationDegree) {
0454       case 1U:
0455         return bins.interpolate1(coords, expDim);
0456       case 3U:
0457         return bins.interpolate3(coords, expDim);
0458       default:
0459         return bins.closest(coords, expDim);
0460     }
0461   }
0462 }  // namespace npstat
0463 
0464 #endif  // NPSTAT_INTERPOLATEHISTOND_HH_