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 }
0468 }
0469
0470 #endif
