Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:02:26 

0001 #ifndef CondFormats_PhysicsToolsObjects_Histogram2D_icc
0002 #define CondFormats_PhysicsToolsObjects_Histogram2D_icc
0003 
0004 #include <cmath>
0005 #include <cstddef>
0006 #include <algorithm>
0007 #include <numeric>
0008 #include <utility>
0009 #include <vector>
0010 
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 
0013 #include "CondFormats/PhysicsToolsObjects/interface/Histogram2D.h"
0014 
0015 namespace PhysicsTools {
0016   namespace Calibration {
0017 
0018     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0019     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D()
0020         : stride(0), limitsX(AxisX_t(), AxisX_t()), limitsY(AxisY_t(), AxisY_t()), total(Value_t()) {
0021       totalValid.store(false, std::memory_order_release);
0022     }
0023 
0024     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0025     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(const Histogram2D &orig)
0026         : stride(orig.stride),
0027           binULimitsX(orig.binULimitsX),
0028           binULimitsY(orig.binULimitsY),
0029           binValues(orig.binValues),
0030           limitsX(orig.limitsX),
0031           limitsY(orig.limitsY),
0032           total(orig.total),
0033           rowTotal(orig.rowTotal),
0034           columnTotal(orig.columnTotal) {
0035       totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0036     }
0037 
0038     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0039     template <typename OValue_t, typename OAxisX_t, typename OAxisY_t>
0040     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(const Histogram2D<OValue_t, OAxisX_t, OAxisY_t> &orig)
0041         : stride(orig.stride),
0042           binULimitsX(orig.binULimitsX.begin(), orig.binULimitsX.end()),
0043           binULimitsY(orig.binULimitsY.begin(), orig.binULimitsY.end()),
0044           binValues(orig.binValues.begin(), orig.binValues.end()),
0045           limitsX(orig.limitsX),
0046           limitsY(orig.limitsY),
0047           total(orig.total),
0048           rowTotal(orig.rowTotal.begin(), orig.rowTotal.end()),
0049           columnTotal(orig.columnTotal.begin(), orig.columnTotal.end()) {
0050       totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0051     }
0052 
0053     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0054     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(const std::vector<AxisX_t> &binULimitsX,
0055                                                         const std::vector<AxisY_t> &binULimitsY)
0056         : stride(binULimitsX.size() + 1),
0057           binULimitsX(binULimitsX),
0058           binULimitsY(binULimitsY),
0059           binValues((binULimitsY.size() + 1) * stride),
0060           limitsX(AxisX_t(), AxisX_t()),
0061           limitsY(AxisY_t(), AxisY_t()),
0062           total(Value_t()),
0063           rowTotal(binULimitsY.size() + 1),
0064           columnTotal(binULimitsX.size() + 1) {
0065       totalValid.store(true, std::memory_order_release);
0066       if (binULimitsX.size() < 2)
0067         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0068                                                    "Fewer than one bin in X requested."
0069                                                 << std::endl;
0070 
0071       limitsX.min = binULimitsX.front();
0072       limitsX.max = binULimitsX.back();
0073 
0074       if (binULimitsY.size() < 2)
0075         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0076                                                    "Fewer than one bin in Y requested."
0077                                                 << std::endl;
0078 
0079       limitsY.min = binULimitsY.front();
0080       limitsY.max = binULimitsY.back();
0081     }
0082 
0083     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0084     template <typename OAxisX_t, typename OAxisY_t>
0085     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(const std::vector<OAxisX_t> &binULimitsX,
0086                                                         const std::vector<OAxisY_t> &binULimitsY)
0087         : stride(binULimitsX.size() + 1),
0088           binULimitsX(binULimitsX.begin(), binULimitsX.end()),
0089           binULimitsY(binULimitsY.begin(), binULimitsY.end()),
0090           binValues((binULimitsY.size() + 1) * stride),
0091           limitsX(AxisX_t(), AxisX_t()),
0092           limitsY(AxisY_t(), AxisY_t()),
0093           total(Value_t()),
0094           rowTotal(binULimitsY.size() + 1),
0095           columnTotal(binULimitsX.size() + 1) {
0096       totalValid.store(true, std::memory_order_release);
0097       if (binULimitsX.size() < 2)
0098         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0099                                                    "Fewer than one bin in X requested."
0100                                                 << std::endl;
0101 
0102       limitsX.min = binULimitsX.front();
0103       limitsX.max = binULimitsX.back();
0104 
0105       if (binULimitsY.size() < 2)
0106         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0107                                                    "Fewer than one bin in Y requested."
0108                                                 << std::endl;
0109 
0110       limitsY.min = binULimitsY.front();
0111       limitsY.max = binULimitsY.back();
0112     }
0113 
0114     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0115     template <typename OAxisX_t, typename OAxisY_t>
0116     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(const std::vector<OAxisX_t> &binULimitsX,
0117                                                         unsigned int nBinsY,
0118                                                         const PhysicsTools::Calibration::Range<OAxisY_t> &rangeY)
0119         : stride(binULimitsX.size() + 1),
0120           binULimitsX(binULimitsX.begin(), binULimitsX.end()),
0121           binValues((nBinsY + 2) * stride),
0122           limitsX(AxisX_t(), AxisX_t()),
0123           limitsY(rangeY),
0124           total(Value_t()),
0125           rowTotal(nBinsY + 2),
0126           columnTotal(binULimitsX.size() + 1) {
0127       totalValid.store(true, std::memory_order_release);
0128       if (binULimitsX.size() < 2)
0129         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0130                                                    "Fewer than one bin in X requested."
0131                                                 << std::endl;
0132 
0133       limitsX.min = binULimitsX.front();
0134       limitsX.max = binULimitsX.back();
0135 
0136       if (!nBinsY)
0137         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0138                                                    "Fewer than one bin in Y requested."
0139                                                 << std::endl;
0140     }
0141 
0142     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0143     template <typename OAxisX_t, typename OAxisY_t>
0144     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(unsigned int nBinsX,
0145                                                         const PhysicsTools::Calibration::Range<OAxisX_t> &rangeX,
0146                                                         const std::vector<OAxisY_t> &binULimitsY)
0147         : stride(nBinsX + 2),
0148           binULimitsY(binULimitsY.begin(), binULimitsY.end()),
0149           binValues((binULimitsY.size() + 1) * stride),
0150           limitsX(rangeX),
0151           limitsY(AxisY_t(), AxisY_t()),
0152           total(Value_t()),
0153           rowTotal(binULimitsY.size() + 1),
0154           columnTotal(nBinsX + 2) {
0155       totalValid.store(true, std::memory_order_release);
0156       if (!nBinsX)
0157         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0158                                                    "Fewer than one bin in X requested."
0159                                                 << std::endl;
0160 
0161       if (binULimitsY.size() < 2)
0162         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0163                                                    "Fewer than one bin in Y requested."
0164                                                 << std::endl;
0165 
0166       limitsY.min = binULimitsY.front();
0167       limitsY.max = binULimitsY.back();
0168     }
0169 
0170     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0171     template <typename OAxisX_t, typename OAxisY_t>
0172     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(unsigned int nBinsX,
0173                                                         const PhysicsTools::Calibration::Range<OAxisX_t> &rangeX,
0174                                                         unsigned int nBinsY,
0175                                                         const PhysicsTools::Calibration::Range<OAxisY_t> &rangeY)
0176         : stride(nBinsX + 2),
0177           binValues((nBinsY + 2) * stride),
0178           limitsX(rangeX),
0179           limitsY(rangeY),
0180           total(Value_t()),
0181           rowTotal(nBinsY + 2),
0182           columnTotal(nBinsX + 2) {
0183       totalValid.store(true, std::memory_order_release);
0184       if (!nBinsX)
0185         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0186                                                    "Fewer than one bin in X requested."
0187                                                 << std::endl;
0188 
0189       if (!nBinsY)
0190         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0191                                                    "Fewer than one bin in Y requested."
0192                                                 << std::endl;
0193     }
0194 
0195     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0196     Histogram2D<Value_t, AxisX_t, AxisY_t>::Histogram2D(
0197         unsigned int nBinsX, AxisX_t minX, AxisX_t maxX, unsigned int nBinsY, AxisY_t minY, AxisY_t maxY)
0198         : stride(nBinsX + 2),
0199           binValues((nBinsY + 2) * stride),
0200           limitsX(minX, maxX),
0201           limitsY(minY, maxY),
0202           total(Value_t()),
0203           rowTotal(nBinsY + 2),
0204           columnTotal(nBinsX + 2) {
0205       totalValid.store(true, std::memory_order_release);
0206       if (!nBinsX)
0207         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0208                                                    "Fewer than one bin in X requested."
0209                                                 << std::endl;
0210 
0211       if (!nBinsY)
0212         throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 2D histogram: "
0213                                                    "Fewer than one bin in Y requested."
0214                                                 << std::endl;
0215     }
0216 
0217     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0218     Histogram2D<Value_t, AxisX_t, AxisY_t>::~Histogram2D() {}
0219 
0220     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0221     Histogram2D<Value_t, AxisX_t, AxisY_t> &Histogram2D<Value_t, AxisX_t, AxisY_t>::operator=(const Histogram2D &orig) {
0222       stride = orig.stride;
0223       binULimitsX = orig.binULimitsX;
0224       binULimitsY = orig.binULimitsY;
0225       binValues = orig.binValues;
0226       limitsX = orig.limitsX;
0227       limitsY = orig.limitsY;
0228       total = orig.total;
0229       totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0230       rowTotal = orig.rowTotal;
0231       columnTotal = orig.columnTotal;
0232       return *this;
0233     }
0234 
0235     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0236     template <typename OValue_t, typename OAxisX_t, typename OAxisY_t>
0237     Histogram2D<Value_t, AxisX_t, AxisY_t> &Histogram2D<Value_t, AxisX_t, AxisY_t>::operator=(
0238         const Histogram2D<OValue_t, OAxisX_t, OAxisY_t> &orig) {
0239       stride = orig.stride;
0240       binULimitsX = std::vector<AxisX_t>(orig.binULimitsX.begin(), orig.binULimitsX.end());
0241       binULimitsY = std::vector<AxisY_t>(orig.binULimitsY.begin(), orig.binULimitsY.end());
0242       binValues = std::vector<Value_t>(orig.binValues.begin(), orig.binValues.end());
0243       limitsX = orig.limitsX;
0244       limitsY = orig.limitsY;
0245       total = orig.total;
0246       totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0247       rowTotal = std::vector<Value_t>(orig.rowTotal.begin(), orig.rowTotal.end());
0248       columnTotal = std::vector<Value_t>(orig.columnTotal.begin(), orig.columnTotal.end());
0249       return *this;
0250     }
0251 
0252     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0253     void Histogram2D<Value_t, AxisX_t, AxisY_t>::reset() {
0254       std::fill(binValues.begin(), binValues.end(), Value_t());
0255       total = Value_t();
0256       totalValid.store(true, std::memory_order_release);
0257       rowTotal.clear();
0258       columnTotal.clear();
0259     }
0260 
0261     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0262     void Histogram2D<Value_t, AxisX_t, AxisY_t>::setBinContent(int bin, Value_t value) {
0263       if (bin < 0 || (unsigned int)bin >= binValues.size())
0264         throw cms::Exception("RangeError") << "Histogram2D bin " << bin
0265                                            << " out of range "
0266                                               "[0, "
0267                                            << (binValues.size() - 1) << "]." << std::endl;
0268 
0269       binValues[bin] = value;
0270       totalValid.store(false, std::memory_order_release);
0271       rowTotal.clear();
0272       columnTotal.clear();
0273     }
0274 
0275     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0276     Value_t Histogram2D<Value_t, AxisX_t, AxisY_t>::normalizedXValue(AxisX_t x, AxisY_t y) const {
0277       int binX = findBinX(x);
0278       int binY = findBinY(y);
0279       return binContent(bin2D(binX, binY)) / normalizationX(binY);
0280     }
0281 
0282     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0283     Value_t Histogram2D<Value_t, AxisX_t, AxisY_t>::normalizedYValue(AxisX_t x, AxisY_t y) const {
0284       int binX = findBinX(x);
0285       int binY = findBinY(y);
0286       return binContent(bin2D(binX, binY)) / normalizationY(binX);
0287     }
0288 
0289     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0290     Value_t Histogram2D<Value_t, AxisX_t, AxisY_t>::normalizedXError(AxisX_t x, AxisY_t y) const {
0291       int binX = findBinX(x);
0292       int binY = findBinY(y);
0293       return std::sqrt(binContent(bin2D(binX, binY))) / normalizationX(binY);
0294     }
0295 
0296     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0297     Value_t Histogram2D<Value_t, AxisX_t, AxisY_t>::normalizedYError(AxisX_t x, AxisY_t y) const {
0298       int binX = findBinX(x);
0299       int binY = findBinY(y);
0300       return std::sqrt(binContent(bin2D(binX, binY))) / normalizationY(binX);
0301     }
0302 
0303     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0304     void Histogram2D<Value_t, AxisX_t, AxisY_t>::fill(AxisX_t x, AxisY_t y, Value_t weight) {
0305       int bin = findBin(x, y);
0306       binValues[bin] += weight;
0307       totalValid.store(false, std::memory_order_release);
0308       rowTotal.clear();
0309       columnTotal.clear();
0310     }
0311 
0312     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0313     void Histogram2D<Value_t, AxisX_t, AxisY_t>::setValues(const std::vector<Value_t> &values) {
0314       if (values.size() != binValues.size())
0315         throw cms::Exception("InvalidVectorSizeError") << "Invalid vector size while setting "
0316                                                           "2D histogram values"
0317                                                        << std::endl;
0318 
0319       binValues = values;
0320       totalValid.store(false, std::memory_order_release);
0321       rowTotal.clear();
0322       columnTotal.clear();
0323     }
0324 
0325     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0326     template <typename OValue_t>
0327     void Histogram2D<Value_t, AxisX_t, AxisY_t>::setValues(const std::vector<OValue_t> &values) {
0328       if (values.size() != binValues.size())
0329         throw cms::Exception("InvalidVectorSizeError") << "Invalid vector size while setting "
0330                                                           "2D histogram values"
0331                                                        << std::endl;
0332 
0333       std::copy(values.begin(), values.end(), binValues.begin());
0334       totalValid.store(false, std::memory_order_release);
0335       rowTotal.clear();
0336       columnTotal.clear();
0337     }
0338 
0339     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0340     typename Histogram2D<Value_t, AxisX_t, AxisY_t>::RangeX Histogram2D<Value_t, AxisX_t, AxisY_t>::binRangeX(
0341         int binX) const {
0342       if (binX < 1 || (unsigned int)binX > stride - 2)
0343         throw cms::Exception("RangeError") << "Histogram2D X bin " << binX << " out of range "
0344                                            << "[1, " << (stride - 2) << "]." << std::endl;
0345 
0346       if (hasEquidistantBinsX()) {
0347         AxisX_t min = (AxisX_t)(binX - 1) / (stride - 2);
0348         AxisX_t max = (AxisX_t)binX / (stride - 2);
0349         min *= limitsX.width();
0350         min += limitsX.min;
0351         max *= limitsX.width();
0352         max += limitsX.min;
0353         return RangeX(min, max);
0354       } else
0355         return RangeX(binULimitsX[binX - 1], binULimitsX[binX]);
0356     }
0357 
0358     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0359     typename Histogram2D<Value_t, AxisX_t, AxisY_t>::RangeY Histogram2D<Value_t, AxisX_t, AxisY_t>::binRangeY(
0360         int binY) const {
0361       unsigned int size = binValues.size() / stride;
0362       if (binY < 1 || (unsigned int)binY > size - 2)
0363         throw cms::Exception("RangeError") << "Histogram2D Y bin " << binY << " out of range "
0364                                            << "[1, " << (size - 2) << "]." << std::endl;
0365 
0366       if (hasEquidistantBinsY()) {
0367         AxisY_t min = (AxisY_t)(binY - 1) / (size - 2);
0368         AxisY_t max = (AxisY_t)binY / (size - 2);
0369         min *= limitsY.width();
0370         min += limitsY.min;
0371         max *= limitsY.width();
0372         max += limitsY.min;
0373         return RangeY(min, max);
0374       } else
0375         return RangeY(binULimitsY[binY - 1], binULimitsY[binY]);
0376     }
0377 
0378     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0379     std::pair<typename Histogram2D<Value_t, AxisX_t, AxisY_t>::RangeX,
0380               typename Histogram2D<Value_t, AxisX_t, AxisY_t>::RangeY>
0381     Histogram2D<Value_t, AxisX_t, AxisY_t>::binRange(int bin) const {
0382       std::size_t size = binValues.size();
0383       if (bin < 1 || (unsigned int)bin > size - 2)
0384         throw cms::Exception("RangeError") << "Histogram2D bin " << bin << " out of range "
0385                                            << "[1, " << (size - 2) << "]." << std::endl;
0386 
0387       return std::make_pair(binRangeX(bin % stride), binRangeY(bin / stride));
0388     }
0389 
0390     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0391     int Histogram2D<Value_t, AxisX_t, AxisY_t>::findBinX(AxisX_t x) const {
0392       if (hasEquidistantBinsX()) {
0393         x -= limitsX.min;
0394         x *= stride - 2;
0395         x /= limitsX.width();
0396         unsigned int iStride = stride;
0397         return clamp(0, (int)(std::floor(x)) + 1, (int)iStride - 1);
0398       } else
0399         return std::upper_bound(binULimitsX.begin(), binULimitsX.end(), x) - binULimitsX.begin();
0400     }
0401 
0402     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0403     int Histogram2D<Value_t, AxisX_t, AxisY_t>::findBinY(AxisY_t y) const {
0404       unsigned int size = binValues.size() / stride;
0405       if (hasEquidistantBinsY()) {
0406         y -= limitsY.min;
0407         y *= size - 2;
0408         y /= limitsY.width();
0409         return clamp(0, (int)(std::floor(y)) + 1, (int)size - 1);
0410       } else
0411         return std::upper_bound(binULimitsY.begin(), binULimitsY.end(), y) - binULimitsY.begin();
0412     }
0413 
0414     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0415     Value_t Histogram2D<Value_t, AxisX_t, AxisY_t>::normalization() const {
0416       if (!totalValid.load(std::memory_order_acquire)) {
0417         total = std::accumulate(binValues.begin() + 1, binValues.end() - 1, Value_t());
0418         totalValid.store(true, std::memory_order_release);
0419       }
0420 
0421       return total;
0422     }
0423 
0424     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0425     Value_t Histogram2D<Value_t, AxisX_t, AxisY_t>::normalizationX(int binY) const {
0426       if (rowTotal.empty()) {
0427         rowTotal.resize(binValues.size() / stride);
0428         typename std::vector<Value_t>::iterator sum = rowTotal.begin();
0429         for (typename std::vector<Value_t>::const_iterator iter = binValues.begin(); iter != binValues.end();
0430              iter += stride)
0431           *sum++ = std::accumulate(iter + 1, iter + (stride - 1), Value_t());
0432       }
0433 
0434       if (binY < 0 || (unsigned int)binY >= binValues.size() / stride)
0435         throw cms::Exception("RangeError") << "Histogram2D bin " << binY
0436                                            << " in Y out of range "
0437                                               "[0, "
0438                                            << (binValues.size() / stride - 1) << "]." << std::endl;
0439 
0440       return rowTotal[binY];
0441     }
0442 
0443     template <typename Value_t, typename AxisX_t, typename AxisY_t>
0444     Value_t Histogram2D<Value_t, AxisX_t, AxisY_t>::normalizationY(int binX) const {
0445       if (columnTotal.empty()) {
0446         columnTotal.resize(stride);
0447         typename std::vector<Value_t>::iterator sum = columnTotal.begin();
0448         for (typename std::vector<Value_t>::const_iterator col = binValues.begin(); col != binValues.begin() + stride;
0449              ++col) {
0450           Value_t value = Value_t();
0451           for (typename std::vector<Value_t>::const_iterator pos = col + stride; pos < (binValues.end() - stride);
0452                pos += stride)
0453             value += *pos;
0454           *sum++ = value;
0455         }
0456       }
0457 
0458       if (binX < 0 || (unsigned int)binX >= stride)
0459         throw cms::Exception("RangeError") << "Histogram2D bin " << binX
0460                                            << " in X out of range "
0461                                               "[0, "
0462                                            << (stride - 1) << "]." << std::endl;
0463 
0464       return columnTotal[binX];
0465     }
0466 
0467   }  // namespace Calibration
0468 }  // namespace PhysicsTools
0469 
0470 #endif  // CondFormats_PhysicsToolsObjects_Histogram2D_icc
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!