PhysicsToolsObjects/interface/Histogram.icc

0001 #ifndef CondFormats_PhysicsToolsObjects_Histogram_icc
0002 #define CondFormats_PhysicsToolsObjects_Histogram_icc
0003
0004 #include <cmath>
0005 #include <cstddef>
0006 #include <algorithm>
0007 #include <numeric>
0008 #include <vector>
0009
0010 #include "FWCore/Utilities/interface/Exception.h"
0011
0012 #include "CondFormats/PhysicsToolsObjects/interface/Histogram.h"
0013
0014 namespace PhysicsTools {
0015   namespace Calibration {
0016
0017     namespace {  // hide when instantiated in -O0 debug mode
0018       template <typename T>
0019       inline const T &clamp(const T &min, const T &val, const T &max) {
0020         if (val < min)
0021           return min;
0022         if (val > max)
0023           return max;
0024         return val;
0025       }
0026     }  // namespace
0027
0028     template <typename Value_t, typename Axis_t>
0029     Histogram<Value_t, Axis_t>::Histogram() : limits(Axis_t(), Axis_t()), total(Value_t()) {
0030       totalValid.store(false, std::memory_order_release);
0031     }
0032
0033     template <typename Value_t, typename Axis_t>
0034     Histogram<Value_t, Axis_t>::Histogram(const Histogram &orig)
0035         : binULimits(orig.binULimits), binValues(orig.binValues), limits(orig.limits), total(orig.total) {
0036       totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0037     }
0038
0039     template <typename Value_t, typename Axis_t>
0040     template <typename OValue_t, typename OAxis_t>
0041     Histogram<Value_t, Axis_t>::Histogram(const Histogram<OValue_t, OAxis_t> &orig)
0042         : binULimits(orig.binULimits.begin(), orig.binULimits.end()),
0043           binValues(orig.binValues.begin(), orig.binValues.end()),
0044           limits(orig.limits),
0045           total(Value_t()) {
0046       totalValid.store(false, std::memory_order_release);
0047     }
0048
0049     template <typename Value_t, typename Axis_t>
0050     Histogram<Value_t, Axis_t>::Histogram(const std::vector<Axis_t> &binULimits)
0051         : binULimits(binULimits), binValues(binULimits.size() + 1), limits(Axis_t(), Axis_t()), total(Value_t()) {
0052       totalValid.store(true, std::memory_order_release);
0053       if (binULimits.size() < 2)
0054         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate histogram: "
0055                                                    "Fewer than one bin requested."
0056                                                 << std::endl;
0057
0058       limits.min = binULimits.front();
0059       limits.max = binULimits.back();
0060     }
0061
0062     template <typename Value_t, typename Axis_t>
0063     template <typename OAxis_t>
0064     Histogram<Value_t, Axis_t>::Histogram(const std::vector<OAxis_t> &binULimits)
0065         : binULimits(binULimits.begin(), binULimits.end()),
0066           binValues(binULimits.size() + 1),
0067           limits(Axis_t(), Axis_t()),
0068           total(Value_t()) {
0069       totalValid.store(true, std::memory_order_release);
0070       if (binULimits.size() < 2)
0071         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate histogram: "
0072                                                    "Fewer than one bin requested."
0073                                                 << std::endl;
0074
0075       limits.min = binULimits.front();
0076       limits.max = binULimits.back();
0077     }
0078
0079     template <typename Value_t, typename Axis_t>
0080     template <typename OAxis_t>
0081     Histogram<Value_t, Axis_t>::Histogram(unsigned int nBins, const PhysicsTools::Calibration::Range<OAxis_t> &range)
0082         : binValues(nBins + 2), limits(range), total(Value_t()) {
0083       totalValid.store(true, std::memory_order_release);
0084       if (!nBins)
0085         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate histogram: "
0086                                                    "Fewer than one bin requested."
0087                                                 << std::endl;
0088     }
0089
0090     template <typename Value_t, typename Axis_t>
0091     Histogram<Value_t, Axis_t>::Histogram(unsigned int nBins, Axis_t min, Axis_t max)
0092         : binValues(nBins + 2), limits(min, max), total(Value_t()) {
0093       totalValid.store(true, std::memory_order_release);
0094       if (!nBins)
0095         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate histogram: "
0096                                                    "Fewer than one bin requested."
0097                                                 << std::endl;
0098     }
0099
0100     template <typename Value_t, typename Axis_t>
0101     Histogram<Value_t, Axis_t>::~Histogram() {}
0102
0103     template <typename Value_t, typename Axis_t>
0104     Histogram<Value_t, Axis_t> &Histogram<Value_t, Axis_t>::operator=(const Histogram &orig) {
0105       binULimits = orig.binULimits;
0106       binValues = orig.binValues;
0107       limits = orig.limits;
0108       total = orig.total;
0109       totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0110       return *this;
0111     }
0112
0113     template <typename Value_t, typename Axis_t>
0114     template <typename OValue_t, typename OAxis_t>
0115     Histogram<Value_t, Axis_t> &Histogram<Value_t, Axis_t>::operator=(const Histogram<OValue_t, OAxis_t> &orig) {
0116       binULimits = std::vector<Axis_t>(orig.binULimits.begin(), orig.binULimits.end());
0117       binValues = std::vector<Value_t>(orig.binValues.begin(), orig.binValues.end());
0118       limits = orig.limits;
0119       total = orig.total;
0120       totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0121       return *this;
0122     }
0123
0124     template <typename Value_t, typename Axis_t>
0125     void Histogram<Value_t, Axis_t>::reset() {
0126       std::fill(binValues.begin(), binValues.end(), Value_t());
0127       total = Value_t();
0128       totalValid.store(true, std::memory_order_release);
0129     }
0130
0131     template <typename Value_t, typename Axis_t>
0132     void Histogram<Value_t, Axis_t>::setBinContent(int bin, Value_t value) {
0133       if (bin < 0 || (unsigned int)bin >= binValues.size())
0134         throw cms::Exception("RangeError") << "Histogram bin " << bin
0135                                            << " out of range "
0136                                               "[0, "
0137                                            << (binValues.size() - 1) << "]." << std::endl;
0138
0139       binValues[bin] = value;
0140       totalValid.store(false, std::memory_order_release);
0141     }
0142
0143     template <typename Value_t, typename Axis_t>
0144     void Histogram<Value_t, Axis_t>::fill(Axis_t x, Value_t weight) {
0145       int bin = findBin(x);
0146       binValues[bin] += weight;
0147       totalValid.store(false, std::memory_order_release);
0148     }
0149
0150     template <typename Value_t, typename Axis_t>
0151     void Histogram<Value_t, Axis_t>::setValues(const std::vector<Value_t> &values) {
0152       if (values.size() != binValues.size())
0153         throw cms::Exception("InvalidVectorSizeError") << "Invalid vector size while setting "
0154                                                           "histogram values"
0155                                                        << std::endl;
0156
0157       binValues = values;
0158       totalValid.store(false, std::memory_order_release);
0159     }
0160
0161     template <typename Value_t, typename Axis_t>
0162     template <typename OValue_t>
0163     void Histogram<Value_t, Axis_t>::setValues(const std::vector<OValue_t> &values) {
0164       if (values.size() != binValues.size())
0165         throw cms::Exception("InvalidVectorSizeError") << "Invalid vector size while setting "
0166                                                           "histogram values"
0167                                                        << std::endl;
0168
0169       std::copy(values.begin(), values.end(), binValues.begin());
0170       totalValid.store(false, std::memory_order_release);
0171     }
0172
0173     template <typename Value_t, typename Axis_t>
0174     typename Histogram<Value_t, Axis_t>::Range Histogram<Value_t, Axis_t>::binRange(int bin) const {
0175       std::size_t size = binValues.size();
0176       if (bin < 1 || (unsigned int)bin > size - 2)
0177         throw cms::Exception("RangeError") << "Histogram bin " << bin << " out of range "
0178                                            << "[1, " << (size - 2) << "]." << std::endl;
0179
0180       if (hasEquidistantBins()) {
0181         Axis_t min = (Axis_t)(bin - 1) / (size - 2);
0182         Axis_t max = (Axis_t)bin / (size - 2);
0183         min *= limits.width();
0184         min += limits.min;
0185         max *= limits.width();
0186         max += limits.min;
0187         return Range(min, max);
0188       } else
0189         return Range(binULimits[bin - 1], binULimits[bin]);
0190     }
0191
0192     template <typename Value_t, typename Axis_t>
0193     int Histogram<Value_t, Axis_t>::findBin(Axis_t x) const {
0194       if (hasEquidistantBins()) {
0195         std::size_t size = binValues.size();
0196         x -= limits.min;
0197         x *= size - 2;
0198         x /= limits.width();
0199         return clamp(0, (int)(std::floor(x)) + 1, (int)size - 1);
0200       } else
0201         return std::upper_bound(binULimits.begin(), binULimits.end(), x) - binULimits.begin();
0202     }
0203
0204     template <typename Value_t, typename Axis_t>
0205     Value_t Histogram<Value_t, Axis_t>::normalization() const {
0206       if (!totalValid.load(std::memory_order_acquire)) {
0207         total = std::accumulate(binValues.begin() + 1, binValues.end() - 1, Value_t());
0208         totalValid.store(true, std::memory_order_release);
0209       }
0210
0211       return total;
0212     }
0213
0214     template <typename Value_t, typename Axis_t>
0215     Value_t Histogram<Value_t, Axis_t>::integral(Axis_t hBound, Axis_t lBound, int mode) const {
0216       if (hBound < lBound)
0217         throw cms::Exception("InvalidIntervalError") << "Upper boundary below lower boundary in "
0218                                                      << "histogram integral." << std::endl;
0219
0220       std::size_t size = binValues.size();
0221       int lBin = clamp(1, findBin(lBound), (int)size - 2);
0222       int hBin = clamp(1, findBin(hBound), (int)size - 2);
0223
0224       double sum = Value_t();
0225       Range lBinRange, hBinRange;
0226
0227       if (hBin > lBin)
0228         sum = std::accumulate(binValues.begin() + (lBin + 1), binValues.begin() + hBin, Value_t());
0229
0230       if (hasEquidistantBins()) {
0231         Axis_t binWidth = limits.width() / (size - 2);
0232         lBinRange = Range((lBin - 1) * binWidth, lBin * binWidth);
0233         hBinRange = Range((hBin - 1) * binWidth, hBin * binWidth);
0234       } else {
0235         lBinRange = Range(binULimits[lBin - 1], binULimits[lBin]);
0236         hBinRange = Range(binULimits[hBin - 1], binULimits[hBin]);
0237       }
0238
0239       switch (mode) {
0240         case 0:
0241           break;
0242         case 1:
0243           lBound = clamp(lBinRange.min, lBound, lBinRange.max);
0244           hBound = clamp(hBinRange.min, hBound, hBinRange.max);
0245           sum += binValues[lBin] * (lBinRange.max - lBound) / lBinRange.width();
0246           sum += binValues[hBin] * (hBound - hBinRange.min) / hBinRange.width();
0247           break;
0248         default:
0249           throw cms::Exception("InvalidMode") << "Invalid mode " << mode
0250                                               << " in "
0251                                                  "Histogram::integral()"
0252                                               << std::endl;
0253       }
0254
0255       return sum;
0256     }
0257
0258   }  // namespace Calibration
0259 }  // namespace PhysicsTools
0260
0261 #endif  // CondFormats_PhysicsToolsObjects_Histogram_icc