InterpolationTables/src/EquidistantSequence.cc

0001 #include <cmath>
0002 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0003
0004 #include "JetMETCorrections/InterpolationTables/interface/EquidistantSequence.h"
0005
0006 namespace npstat {
0007   EquidistantInLinearSpace::EquidistantInLinearSpace(const double minScale,
0008                                                      const double maxScale,
0009                                                      const unsigned nScales)
0010       : std::vector<double>() {
0011     switch (nScales) {
0012       case 0:
0013         break;
0014
0015       case 1: {
0016         this->reserve(nScales);
0017         const double sc = (minScale == maxScale ? minScale : (minScale + maxScale) / 2.0);
0018         push_back(sc);
0019       } break;
0020
0021       default: {
0022         this->reserve(nScales);
0023         const double step = (maxScale - minScale) / (nScales - 1);
0024         push_back(minScale);
0025         for (unsigned i = 1; i < nScales - 1; ++i)
0026           push_back(minScale + i * step);
0027         push_back(maxScale);
0028       } break;
0029     }
0030   }
0031
0032   EquidistantInLogSpace::EquidistantInLogSpace(const double minScale, const double maxScale, const unsigned nScales)
0033       : std::vector<double>() {
0034     if (nScales)
0035       if (!(minScale > 0.0 && maxScale > 0.0))
0036         throw npstat::NpstatInvalidArgument(
0037             "In npstat::EquidistantInLogSpace constructor: "
0038             "minimum and maximum scales must be positive");
0039     switch (nScales) {
0040       case 0:
0041         break;
0042
0043       case 1: {
0044         this->reserve(nScales);
0045         const double sc = (minScale == maxScale ? minScale : sqrt(minScale * maxScale));
0046         push_back(sc);
0047       } break;
0048
0049       default: {
0050         this->reserve(nScales);
0051         const double logmax = log(maxScale);
0052         const double logmin = log(minScale);
0053         const double logstep = (logmax - logmin) / (nScales - 1);
0054         push_back(minScale);
0055         for (unsigned i = 1; i < nScales - 1; ++i)
0056           push_back(exp(logmin + i * logstep));
0057         push_back(maxScale);
0058       } break;
0059     }
0060   }
0061 }  // namespace npstat