00001 #ifndef NPSTAT_INTERPOLATEHISTOND_HH_ 00002 #define NPSTAT_INTERPOLATEHISTOND_HH_ 00003 00024 #include "JetMETCorrections/InterpolationTables/interface/HistoND.h" 00025 00026 namespace npstat { 00034 template <typename Float, class Axis> 00035 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00036 const double *coords, unsigned coordsDim, 00037 unsigned interpolationDegree); 00039 00043 template <typename Float, class Axis> 00044 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00045 double x0, unsigned interpolationDegree); 00046 00047 template <typename Float, class Axis> 00048 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00049 double x0, double x1, 00050 unsigned interpolationDegree); 00051 00052 template <typename Float, class Axis> 00053 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00054 double x0, double x1, double x2, 00055 unsigned interpolationDegree); 00056 00057 template <typename Float, class Axis> 00058 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00059 double x0, double x1, double x2, double x3, 00060 unsigned interpolationDegree); 00061 00062 template <typename Float, class Axis> 00063 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00064 double x0, double x1, double x2, double x3, 00065 double x4, unsigned interpolationDegree); 00066 00067 template <typename Float, class Axis> 00068 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00069 double x0, double x1, double x2, double x3, 00070 double x4, double x5, 00071 unsigned interpolationDegree); 00072 00073 template <typename Float, class Axis> 00074 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00075 double x0, double x1, double x2, double x3, 00076 double x4, double x5, double x6, 00077 unsigned interpolationDegree); 00078 00079 template <typename Float, class Axis> 00080 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00081 double x0, double x1, double x2, double x3, 00082 double x4, double x5, double x6, double x7, 00083 unsigned interpolationDegree); 00084 00085 template <typename Float, class Axis> 00086 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00087 double x0, double x1, double x2, double x3, 00088 double x4, double x5, double x6, double x7, 00089 double x8, unsigned interpolationDegree); 00090 00091 template <typename Float, class Axis> 00092 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00093 double x0, double x1, double x2, double x3, 00094 double x4, double x5, double x6, double x7, 00095 double x8, double x9, 00096 unsigned interpolationDegree); 00098 } 00099 00100 #include <cassert> 00101 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h" 00102 00103 namespace npstat { 00104 namespace Private { 00105 template <typename Float, class Axis> 00106 void iHND_checkArgs(const HistoND<Float,Axis>& histo, 00107 const unsigned xDim, 00108 const unsigned interpolationDegree) 00109 { 00110 if (xDim != histo.dim()) throw npstat::NpstatInvalidArgument( 00111 "In npstat::interpolateHistoND: incompatible " 00112 "dimensionality of input coordinates"); 00113 if (xDim == 0U) throw npstat::NpstatInvalidArgument( 00114 "In npstat::interpolateHistoND: can not interpolate " 00115 "zero-dimensional histograms"); 00116 if (!(interpolationDegree == 0U || 00117 interpolationDegree == 1U || 00118 interpolationDegree == 3U)) throw npstat::NpstatInvalidArgument( 00119 "In npstat::interpolateHistoND: " 00120 "unsupported interpolation degree"); 00121 if (interpolationDegree == 3U && !histo.isUniformlyBinned()) 00122 throw npstat::NpstatInvalidArgument( 00123 "In npstat::interpolateHistoND: unsupported " 00124 "interpolation degree for non-uniform binning"); 00125 } 00126 } 00127 00128 template <typename Float, class Axis> 00129 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00130 const double *x, const unsigned xDim, 00131 const unsigned interpolationDegree) 00132 { 00133 Private::iHND_checkArgs(histo, xDim, interpolationDegree); 00134 assert(x); 00135 const Axis* ax = &histo.axes()[0]; 00136 double coords[CHAR_BIT*sizeof(unsigned long)]; 00137 for (unsigned i=0; i<xDim; ++i) 00138 coords[i] = ax[i].fltBinNumber(x[i], false); 00139 const ArrayND<Float>& bins(histo.binContents()); 00140 switch (interpolationDegree) 00141 { 00142 case 1U: 00143 return bins.interpolate1(coords, xDim); 00144 case 3U: 00145 return bins.interpolate3(coords, xDim); 00146 default: 00147 return bins.closest(coords, xDim); 00148 } 00149 } 00150 00151 template <typename Float, class Axis> 00152 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00153 const double x0, 00154 const unsigned interpolationDegree) 00155 { 00156 const unsigned expDim = 1U; 00157 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00158 const double coords = histo.axis(0).fltBinNumber(x0, false); 00159 const ArrayND<Float>& bins(histo.binContents()); 00160 switch (interpolationDegree) 00161 { 00162 case 1U: 00163 return bins.interpolate1(&coords, expDim); 00164 case 3U: 00165 return bins.interpolate3(&coords, expDim); 00166 default: 00167 return bins.closest(&coords, expDim); 00168 } 00169 } 00170 00171 template <typename Float, class Axis> 00172 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00173 const double x0, 00174 const double x1, 00175 const unsigned interpolationDegree) 00176 { 00177 const unsigned expDim = 2U; 00178 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00179 const Axis* ax = &histo.axes()[0]; 00180 double coords[expDim]; 00181 coords[0] = ax[0].fltBinNumber(x0, false); 00182 coords[1] = ax[1].fltBinNumber(x1, false); 00183 const ArrayND<Float>& bins(histo.binContents()); 00184 switch (interpolationDegree) 00185 { 00186 case 1U: 00187 return bins.interpolate1(coords, expDim); 00188 case 3U: 00189 return bins.interpolate3(coords, expDim); 00190 default: 00191 return bins.closest(coords, expDim); 00192 } 00193 } 00194 00195 template <typename Float, class Axis> 00196 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00197 const double x0, 00198 const double x1, 00199 const double x2, 00200 const unsigned interpolationDegree) 00201 { 00202 const unsigned expDim = 3U; 00203 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00204 const Axis* ax = &histo.axes()[0]; 00205 double coords[expDim]; 00206 coords[0] = ax[0].fltBinNumber(x0, false); 00207 coords[1] = ax[1].fltBinNumber(x1, false); 00208 coords[2] = ax[2].fltBinNumber(x2, false); 00209 const ArrayND<Float>& bins(histo.binContents()); 00210 switch (interpolationDegree) 00211 { 00212 case 1U: 00213 return bins.interpolate1(coords, expDim); 00214 case 3U: 00215 return bins.interpolate3(coords, expDim); 00216 default: 00217 return bins.closest(coords, expDim); 00218 } 00219 } 00220 00221 template <typename Float, class Axis> 00222 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00223 const double x0, 00224 const double x1, 00225 const double x2, 00226 const double x3, 00227 const unsigned interpolationDegree) 00228 { 00229 const unsigned expDim = 4U; 00230 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00231 const Axis* ax = &histo.axes()[0]; 00232 double coords[expDim]; 00233 coords[0] = ax[0].fltBinNumber(x0, false); 00234 coords[1] = ax[1].fltBinNumber(x1, false); 00235 coords[2] = ax[2].fltBinNumber(x2, false); 00236 coords[3] = ax[3].fltBinNumber(x3, false); 00237 const ArrayND<Float>& bins(histo.binContents()); 00238 switch (interpolationDegree) 00239 { 00240 case 1U: 00241 return bins.interpolate1(coords, expDim); 00242 case 3U: 00243 return bins.interpolate3(coords, expDim); 00244 default: 00245 return bins.closest(coords, expDim); 00246 } 00247 } 00248 00249 template <typename Float, class Axis> 00250 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00251 const double x0, 00252 const double x1, 00253 const double x2, 00254 const double x3, 00255 const double x4, 00256 const unsigned interpolationDegree) 00257 { 00258 const unsigned expDim = 5U; 00259 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00260 const Axis* ax = &histo.axes()[0]; 00261 double coords[expDim]; 00262 coords[0] = ax[0].fltBinNumber(x0, false); 00263 coords[1] = ax[1].fltBinNumber(x1, false); 00264 coords[2] = ax[2].fltBinNumber(x2, false); 00265 coords[3] = ax[3].fltBinNumber(x3, false); 00266 coords[4] = ax[4].fltBinNumber(x4, false); 00267 const ArrayND<Float>& bins(histo.binContents()); 00268 switch (interpolationDegree) 00269 { 00270 case 1U: 00271 return bins.interpolate1(coords, expDim); 00272 case 3U: 00273 return bins.interpolate3(coords, expDim); 00274 default: 00275 return bins.closest(coords, expDim); 00276 } 00277 } 00278 00279 template <typename Float, class Axis> 00280 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00281 const double x0, 00282 const double x1, 00283 const double x2, 00284 const double x3, 00285 const double x4, 00286 const double x5, 00287 const unsigned interpolationDegree) 00288 { 00289 const unsigned expDim = 6U; 00290 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00291 const Axis* ax = &histo.axes()[0]; 00292 double coords[expDim]; 00293 coords[0] = ax[0].fltBinNumber(x0, false); 00294 coords[1] = ax[1].fltBinNumber(x1, false); 00295 coords[2] = ax[2].fltBinNumber(x2, false); 00296 coords[3] = ax[3].fltBinNumber(x3, false); 00297 coords[4] = ax[4].fltBinNumber(x4, false); 00298 coords[5] = ax[5].fltBinNumber(x5, false); 00299 const ArrayND<Float>& bins(histo.binContents()); 00300 switch (interpolationDegree) 00301 { 00302 case 1U: 00303 return bins.interpolate1(coords, expDim); 00304 case 3U: 00305 return bins.interpolate3(coords, expDim); 00306 default: 00307 return bins.closest(coords, expDim); 00308 } 00309 } 00310 00311 template <typename Float, class Axis> 00312 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00313 const double x0, 00314 const double x1, 00315 const double x2, 00316 const double x3, 00317 const double x4, 00318 const double x5, 00319 const double x6, 00320 const unsigned interpolationDegree) 00321 { 00322 const unsigned expDim = 7U; 00323 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00324 const Axis* ax = &histo.axes()[0]; 00325 double coords[expDim]; 00326 coords[0] = ax[0].fltBinNumber(x0, false); 00327 coords[1] = ax[1].fltBinNumber(x1, false); 00328 coords[2] = ax[2].fltBinNumber(x2, false); 00329 coords[3] = ax[3].fltBinNumber(x3, false); 00330 coords[4] = ax[4].fltBinNumber(x4, false); 00331 coords[5] = ax[5].fltBinNumber(x5, false); 00332 coords[6] = ax[6].fltBinNumber(x6, false); 00333 const ArrayND<Float>& bins(histo.binContents()); 00334 switch (interpolationDegree) 00335 { 00336 case 1U: 00337 return bins.interpolate1(coords, expDim); 00338 case 3U: 00339 return bins.interpolate3(coords, expDim); 00340 default: 00341 return bins.closest(coords, expDim); 00342 } 00343 } 00344 00345 template <typename Float, class Axis> 00346 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00347 const double x0, 00348 const double x1, 00349 const double x2, 00350 const double x3, 00351 const double x4, 00352 const double x5, 00353 const double x6, 00354 const double x7, 00355 const unsigned interpolationDegree) 00356 { 00357 const unsigned expDim = 8U; 00358 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00359 const Axis* ax = &histo.axes()[0]; 00360 double coords[expDim]; 00361 coords[0] = ax[0].fltBinNumber(x0, false); 00362 coords[1] = ax[1].fltBinNumber(x1, false); 00363 coords[2] = ax[2].fltBinNumber(x2, false); 00364 coords[3] = ax[3].fltBinNumber(x3, false); 00365 coords[4] = ax[4].fltBinNumber(x4, false); 00366 coords[5] = ax[5].fltBinNumber(x5, false); 00367 coords[6] = ax[6].fltBinNumber(x6, false); 00368 coords[7] = ax[7].fltBinNumber(x7, false); 00369 const ArrayND<Float>& bins(histo.binContents()); 00370 switch (interpolationDegree) 00371 { 00372 case 1U: 00373 return bins.interpolate1(coords, expDim); 00374 case 3U: 00375 return bins.interpolate3(coords, expDim); 00376 default: 00377 return bins.closest(coords, expDim); 00378 } 00379 } 00380 00381 template <typename Float, class Axis> 00382 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00383 const double x0, 00384 const double x1, 00385 const double x2, 00386 const double x3, 00387 const double x4, 00388 const double x5, 00389 const double x6, 00390 const double x7, 00391 const double x8, 00392 const unsigned interpolationDegree) 00393 { 00394 const unsigned expDim = 9U; 00395 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00396 const Axis* ax = &histo.axes()[0]; 00397 double coords[expDim]; 00398 coords[0] = ax[0].fltBinNumber(x0, false); 00399 coords[1] = ax[1].fltBinNumber(x1, false); 00400 coords[2] = ax[2].fltBinNumber(x2, false); 00401 coords[3] = ax[3].fltBinNumber(x3, false); 00402 coords[4] = ax[4].fltBinNumber(x4, false); 00403 coords[5] = ax[5].fltBinNumber(x5, false); 00404 coords[6] = ax[6].fltBinNumber(x6, false); 00405 coords[7] = ax[7].fltBinNumber(x7, false); 00406 coords[8] = ax[8].fltBinNumber(x8, false); 00407 const ArrayND<Float>& bins(histo.binContents()); 00408 switch (interpolationDegree) 00409 { 00410 case 1U: 00411 return bins.interpolate1(coords, expDim); 00412 case 3U: 00413 return bins.interpolate3(coords, expDim); 00414 default: 00415 return bins.closest(coords, expDim); 00416 } 00417 } 00418 00419 template <typename Float, class Axis> 00420 Float interpolateHistoND(const HistoND<Float,Axis>& histo, 00421 const double x0, 00422 const double x1, 00423 const double x2, 00424 const double x3, 00425 const double x4, 00426 const double x5, 00427 const double x6, 00428 const double x7, 00429 const double x8, 00430 const double x9, 00431 const unsigned interpolationDegree) 00432 { 00433 const unsigned expDim = 10U; 00434 Private::iHND_checkArgs(histo, expDim, interpolationDegree); 00435 const Axis* ax = &histo.axes()[0]; 00436 double coords[expDim]; 00437 coords[0] = ax[0].fltBinNumber(x0, false); 00438 coords[1] = ax[1].fltBinNumber(x1, false); 00439 coords[2] = ax[2].fltBinNumber(x2, false); 00440 coords[3] = ax[3].fltBinNumber(x3, false); 00441 coords[4] = ax[4].fltBinNumber(x4, false); 00442 coords[5] = ax[5].fltBinNumber(x5, false); 00443 coords[6] = ax[6].fltBinNumber(x6, false); 00444 coords[7] = ax[7].fltBinNumber(x7, false); 00445 coords[8] = ax[8].fltBinNumber(x8, false); 00446 coords[9] = ax[9].fltBinNumber(x9, false); 00447 const ArrayND<Float>& bins(histo.binContents()); 00448 switch (interpolationDegree) 00449 { 00450 case 1U: 00451 return bins.interpolate1(coords, expDim); 00452 case 3U: 00453 return bins.interpolate3(coords, expDim); 00454 default: 00455 return bins.closest(coords, expDim); 00456 } 00457 } 00458 } 00459 00460 00461 #endif // NPSTAT_INTERPOLATEHISTOND_HH_ 00462