File indexing completed on 2024-04-06 12:02:26
0001 #ifndef CondFormats_PhysicsToolsObjects_Histogram3D_icc
0002 #define CondFormats_PhysicsToolsObjects_Histogram3D_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/Histogram3D.h"
0014
0015 namespace PhysicsTools {
0016 namespace Calibration {
0017
0018 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0019 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::Histogram3D()
0020 : strideX(0),
0021 strideY(0),
0022 limitsX(AxisX_t(), AxisX_t()),
0023 limitsY(AxisY_t(), AxisY_t()),
0024 limitsZ(AxisZ_t(), AxisZ_t()),
0025 total(Value_t()) {
0026 totalValid.store(false, std::memory_order_release);
0027 }
0028
0029 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0030 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::Histogram3D(const Histogram3D &orig)
0031 : strideX(orig.strideX),
0032 strideY(orig.strideY),
0033 binULimitsX(orig.binULimitsX),
0034 binULimitsY(orig.binULimitsY),
0035 binULimitsZ(orig.binULimitsZ),
0036 binValues(orig.binValues),
0037 limitsX(orig.limitsX),
0038 limitsY(orig.limitsY),
0039 limitsZ(orig.limitsZ),
0040 total(orig.total),
0041 rowTotal(orig.rowTotal),
0042 columnTotal(orig.columnTotal) {
0043 totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0044 }
0045
0046 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0047 template <typename OValue_t, typename OAxisX_t, typename OAxisY_t, typename OAxisZ_t>
0048 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::Histogram3D(
0049 const Histogram3D<OValue_t, OAxisX_t, OAxisY_t, OAxisZ_t> &orig)
0050 : strideX(orig.strideX),
0051 strideY(orig.strideY),
0052 binULimitsX(orig.binULimitsX.begin(), orig.binULimitsX.end()),
0053 binULimitsY(orig.binULimitsY.begin(), orig.binULimitsY.end()),
0054 binULimitsZ(orig.binULimitsZ.begin(), orig.binULimitsZ.end()),
0055 binValues(orig.binValues.begin(), orig.binValues.end()),
0056 limitsX(orig.limitsX),
0057 limitsY(orig.limitsY),
0058 limitsZ(orig.limitsZ),
0059 total(orig.total),
0060 rowTotal(orig.rowTotal.begin(), orig.rowTotal.end()),
0061 columnTotal(orig.columnTotal.begin(), orig.columnTotal.end()) {
0062 totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0063 }
0064
0065 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0066 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::Histogram3D(const std::vector<AxisX_t> &binULimitsX,
0067 const std::vector<AxisY_t> &binULimitsY,
0068 const std::vector<AxisZ_t> &binULimitsZ)
0069 : strideX(binULimitsX.size() + 1),
0070 strideY(binULimitsY.size() + 1),
0071 binULimitsX(binULimitsX),
0072 binULimitsY(binULimitsY),
0073 binULimitsZ(binULimitsZ),
0074 binValues((binULimitsZ.size() + 1) * strideX * strideY),
0075 limitsX(AxisX_t(), AxisX_t()),
0076 limitsY(AxisY_t(), AxisY_t()),
0077 limitsZ(AxisZ_t(), AxisZ_t()),
0078 total(Value_t()),
0079 sliceTotal(binULimitsZ.size() + 1),
0080 rowTotal(binULimitsY.size() + 1),
0081 columnTotal(binULimitsX.size() + 1) {
0082 totalValid.store(true, std::memory_order_release);
0083 if (binULimitsX.size() < 2)
0084 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0085 "Fewer than one bin in X requested."
0086 << std::endl;
0087
0088 limitsX.min = binULimitsX.front();
0089 limitsX.max = binULimitsX.back();
0090
0091 if (binULimitsY.size() < 2)
0092 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0093 "Fewer than one bin in Y requested."
0094 << std::endl;
0095
0096 limitsY.min = binULimitsY.front();
0097 limitsY.max = binULimitsY.back();
0098
0099 if (binULimitsZ.size() < 2)
0100 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0101 "Fewer than one bin in 3 requested."
0102 << std::endl;
0103
0104 limitsZ.min = binULimitsZ.front();
0105 limitsZ.max = binULimitsZ.back();
0106 }
0107
0108 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0109 template <typename OAxisX_t, typename OAxisY_t, typename OAxisZ_t>
0110 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::Histogram3D(const std::vector<OAxisX_t> &binULimitsX,
0111 const std::vector<OAxisY_t> &binULimitsY,
0112 const std::vector<OAxisZ_t> &binULimitsZ)
0113 : strideX(binULimitsX.size() + 1),
0114 strideY(binULimitsY.size() + 1),
0115 binULimitsX(binULimitsX.begin(), binULimitsX.end()),
0116 binULimitsY(binULimitsY.begin(), binULimitsY.end()),
0117 binULimitsZ(binULimitsZ.begin(), binULimitsZ.end()),
0118 binValues((binULimitsZ.size() + 1) * strideX * strideY),
0119 limitsX(AxisX_t(), AxisX_t()),
0120 limitsY(AxisY_t(), AxisY_t()),
0121 limitsZ(AxisZ_t(), AxisZ_t()),
0122 total(Value_t()),
0123 sliceTotal(binULimitsY.size() + 1),
0124 rowTotal(binULimitsY.size() + 1),
0125 columnTotal(binULimitsX.size() + 1) {
0126 totalValid.store(true, std::memory_order_release);
0127 if (binULimitsX.size() < 2)
0128 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0129 "Fewer than one bin in X requested."
0130 << std::endl;
0131
0132 limitsX.min = binULimitsX.front();
0133 limitsX.max = binULimitsX.back();
0134
0135 if (binULimitsY.size() < 2)
0136 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0137 "Fewer than one bin in Y requested."
0138 << std::endl;
0139
0140 limitsY.min = binULimitsY.front();
0141 limitsY.max = binULimitsY.back();
0142
0143 if (binULimitsZ.size() < 2)
0144 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0145 "Fewer than one bin in Z requested."
0146 << std::endl;
0147
0148 limitsZ.min = binULimitsZ.front();
0149 limitsZ.max = binULimitsZ.back();
0150 }
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0232 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::Histogram3D(unsigned int nBinsX,
0233 AxisX_t minX,
0234 AxisX_t maxX,
0235 unsigned int nBinsY,
0236 AxisY_t minY,
0237 AxisY_t maxY,
0238 unsigned int nBinsZ,
0239 AxisZ_t minZ,
0240 AxisZ_t maxZ)
0241 : strideX(nBinsX + 2),
0242 strideY(nBinsY + 2),
0243 binValues((nBinsZ + 2) * strideX * strideY),
0244 limitsX(minX, maxX),
0245 limitsY(minY, maxY),
0246 limitsZ(minZ, maxZ),
0247 total(Value_t()),
0248 sliceTotal(nBinsZ + 2),
0249 rowTotal(nBinsY + 2),
0250 columnTotal(nBinsX + 2) {
0251 totalValid.store(true, std::memory_order_release);
0252 if (!nBinsX)
0253 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0254 "Fewer than one bin in X requested."
0255 << std::endl;
0256
0257 if (!nBinsY)
0258 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0259 "Fewer than one bin in Y requested."
0260 << std::endl;
0261
0262 if (!nBinsZ)
0263 throw cms::Exception("TooFewBinsError") << "Trying to generate degenerate 3D histogram: "
0264 "Fewer than one bin in Z requested."
0265 << std::endl;
0266 }
0267
0268 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0269 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::~Histogram3D() {}
0270
0271 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0272 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t> &Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::operator=(
0273 const Histogram3D &orig) {
0274 strideX = orig.strideX;
0275 strideY = orig.strideY;
0276 binULimitsX = orig.binULimitsX;
0277 binULimitsY = orig.binULimitsY;
0278 binULimitsZ = orig.binULimitsZ;
0279 binValues = orig.binValues;
0280 limitsX = orig.limitsX;
0281 limitsY = orig.limitsY;
0282 limitsZ = orig.limitsZ;
0283 total = orig.total;
0284 totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0285 rowTotal = orig.rowTotal;
0286 columnTotal = orig.columnTotal;
0287 sliceTotal = orig.sliceTotal;
0288 return *this;
0289 }
0290
0291 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0292 template <typename OValue_t, typename OAxisX_t, typename OAxisY_t, typename OAxisZ_t>
0293 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t> &Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::operator=(
0294 const Histogram3D<OValue_t, OAxisX_t, OAxisY_t, OAxisZ_t> &orig) {
0295 strideX = orig.strideX;
0296 strideY = orig.strideY;
0297 binULimitsX = std::vector<AxisX_t>(orig.binULimitsX.begin(), orig.binULimitsX.end());
0298 binULimitsY = std::vector<AxisY_t>(orig.binULimitsY.begin(), orig.binULimitsY.end());
0299 binULimitsZ = std::vector<AxisZ_t>(orig.binULimitsZ.begin(), orig.binULimitsZ.end());
0300 binValues = std::vector<Value_t>(orig.binValues.begin(), orig.binValues.end());
0301 limitsX = orig.limitsX;
0302 limitsY = orig.limitsY;
0303 limitsZ = orig.limitsZ;
0304 total = orig.total;
0305 totalValid.store(orig.totalValid.load(std::memory_order_acquire), std::memory_order_release);
0306 rowTotal = std::vector<Value_t>(orig.rowTotal.begin(), orig.rowTotal.end());
0307 columnTotal = std::vector<Value_t>(orig.columnTotal.begin(), orig.columnTotal.end());
0308 sliceTotal = std::vector<Value_t>(orig.sliceTotal.begin(), orig.sliceTotal.end());
0309
0310 return *this;
0311 }
0312
0313 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0314 void Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::reset() {
0315 std::fill(binValues.begin(), binValues.end(), Value_t());
0316 total = Value_t();
0317 totalValid.store(true, std::memory_order_release);
0318 sliceTotal.clear();
0319 rowTotal.clear();
0320 columnTotal.clear();
0321 }
0322
0323 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0324 void Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::setBinContent(int bin, Value_t value) {
0325 if (bin < 0 || (unsigned int)bin >= binValues.size())
0326 throw cms::Exception("RangeError") << "Histogram3D bin " << bin
0327 << " out of range "
0328 "[0, "
0329 << (binValues.size() - 1) << "]." << std::endl;
0330
0331 binValues[bin] = value;
0332 totalValid.store(false, std::memory_order_release);
0333 rowTotal.clear();
0334 columnTotal.clear();
0335 sliceTotal.clear();
0336 }
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0379 void Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::fill(AxisX_t x, AxisY_t y, AxisZ_t z, Value_t weight) {
0380 int bin = findBin(x, y, z);
0381 binValues[bin] += weight;
0382 totalValid.store(false, std::memory_order_release);
0383 rowTotal.clear();
0384 columnTotal.clear();
0385 sliceTotal.clear();
0386 }
0387
0388 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0389 void Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::setValues(const std::vector<Value_t> &values) {
0390 if (values.size() != binValues.size())
0391 throw cms::Exception("InvalidVectorSizeError") << "Invalid vector size while setting "
0392 "3D histogram values"
0393 << std::endl;
0394
0395 binValues = values;
0396 totalValid.store(false, std::memory_order_release);
0397 rowTotal.clear();
0398 columnTotal.clear();
0399 sliceTotal.clear();
0400 }
0401
0402 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0403 template <typename OValue_t>
0404 void Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::setValues(const std::vector<OValue_t> &values) {
0405 if (values.size() != binValues.size())
0406 throw cms::Exception("InvalidVectorSizeError") << "Invalid vector size while setting "
0407 "3D histogram values"
0408 << std::endl;
0409
0410 std::copy(values.begin(), values.end(), binValues.begin());
0411 totalValid.store(false, std::memory_order_release);
0412 rowTotal.clear();
0413 columnTotal.clear();
0414 sliceTotal.clear();
0415 }
0416
0417 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0418 typename Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::RangeX
0419 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::binRangeX(int binX) const {
0420 if (binX < 1 || (unsigned int)binX > strideX - 2)
0421 throw cms::Exception("RangeError") << "Histogram3D X bin " << binX << " out of range "
0422 << "[1, " << (strideX - 2) << "]." << std::endl;
0423
0424 if (hasEquidistantBinsX()) {
0425 AxisX_t min = (AxisX_t)(binX - 1) / (strideX - 2);
0426 AxisX_t max = (AxisX_t)binX / (strideX - 2);
0427 min *= limitsX.width();
0428 min += limitsX.min;
0429 max *= limitsX.width();
0430 max += limitsX.min;
0431 return RangeX(min, max);
0432 } else
0433 return RangeX(binULimitsX[binX - 1], binULimitsX[binX]);
0434 }
0435
0436 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0437 typename Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::RangeY
0438 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::binRangeY(int binY) const {
0439 if (binY < 1 || (unsigned int)binY > strideY - 2)
0440 throw cms::Exception("RangeError") << "Histogram3D Y bin " << binY << " out of range "
0441 << "[1, " << (strideY - 2) << "]." << std::endl;
0442
0443 if (hasEquidistantBinsY()) {
0444 AxisY_t min = (AxisY_t)(binY - 1) / (strideY - 2);
0445 AxisY_t max = (AxisY_t)binY / (strideY - 2);
0446 min *= limitsY.width();
0447 min += limitsY.min;
0448 max *= limitsY.width();
0449 max += limitsY.min;
0450 return RangeY(min, max);
0451 } else
0452 return RangeY(binULimitsY[binY - 1], binULimitsY[binY]);
0453 }
0454
0455 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0456 typename Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::RangeZ
0457 Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::binRangeZ(int binZ) const {
0458 unsigned int size = binValues.size() / (strideX * strideY);
0459 if (binZ < 1 || (unsigned int)binZ > size - 2)
0460 throw cms::Exception("RangeError") << "Histogram3D Z bin " << binZ << " out of range "
0461 << "[1, " << (size - 2) << "]." << std::endl;
0462
0463 if (hasEquidistantBinsZ()) {
0464 AxisZ_t min = (AxisY_t)(binZ - 1) / (size - 2);
0465 AxisZ_t max = (AxisY_t)binZ / (size - 2);
0466 min *= limitsZ.width();
0467 min += limitsZ.min;
0468 max *= limitsZ.width();
0469 max += limitsZ.min;
0470 return RangeZ(min, max);
0471 } else
0472 return RangeZ(binULimitsZ[binZ - 1], binULimitsZ[binZ]);
0473 }
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0495 int Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::findBinX(AxisX_t x) const {
0496 if (hasEquidistantBinsX()) {
0497 x -= limitsX.min;
0498 x *= strideX - 2;
0499 x /= limitsX.width();
0500 unsigned int iStrideX = strideX;
0501 return clamp(0, (int)(std::floor(x)) + 1, (int)iStrideX - 1);
0502 } else
0503 return std::upper_bound(binULimitsX.begin(), binULimitsX.end(), x) - binULimitsX.begin();
0504 }
0505
0506 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0507 int Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::findBinY(AxisY_t y) const {
0508 if (hasEquidistantBinsY()) {
0509 y -= limitsY.min;
0510 y *= strideY - 2;
0511 y /= limitsY.width();
0512 unsigned int iStrideY = strideY;
0513 return clamp(0, (int)(std::floor(y)) + 1, (int)iStrideY - 1);
0514 } else
0515 return std::upper_bound(binULimitsY.begin(), binULimitsY.end(), y) - binULimitsY.begin();
0516 }
0517
0518 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0519 int Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::findBinZ(AxisZ_t z) const {
0520 unsigned int size = binValues.size() / (strideX * strideY);
0521 if (hasEquidistantBinsZ()) {
0522 z -= limitsZ.min;
0523 z *= size - 2;
0524 z /= limitsZ.width();
0525 return clamp(0, (int)(std::floor(z)) + 1, (int)size - 1);
0526 } else
0527 return std::upper_bound(binULimitsZ.begin(), binULimitsZ.end(), z) - binULimitsZ.begin();
0528 }
0529
0530 template <typename Value_t, typename AxisX_t, typename AxisY_t, typename AxisZ_t>
0531 Value_t Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::normalization() const {
0532 if (!totalValid.load(std::memory_order_acquire)) {
0533 total = std::accumulate(binValues.begin() + 1, binValues.end() - 1, Value_t());
0534 totalValid.store(true, std::memory_order_release);
0535 }
0536
0537 return total;
0538 }
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591 }
0592 }
0593
0594 #endif