File indexing completed on 2023-03-17 11:10:41
0001 #ifndef NPSTAT_INTERPOLATEHISTOND_HH_
0002 #define NPSTAT_INTERPOLATEHISTOND_HH_
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 #include "JetMETCorrections/InterpolationTables/interface/HistoND.h"
0025
0026 namespace npstat {
0027
0028
0029
0030
0031
0032
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
0042
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 }
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 }
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 }
0463
0464 #endif