File indexing completed on 2024-04-06 12:02:15
0001 #include <cmath>
0002 #include <climits>
0003
0004 #include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h"
0005
0006 HBHENegativeEFilter::HBHENegativeEFilter(const std::vector<PiecewiseScalingPolynomial>& a1vec,
0007 const std::vector<PiecewiseScalingPolynomial>& a2vec,
0008 const std::vector<uint32_t>& iEtaLimits,
0009 const std::vector<std::pair<double, double> >& cut,
0010 const double minCharge,
0011 const unsigned firstTimeSlice,
0012 const unsigned lastTimeSlice)
0013 : a1v_(a1vec),
0014 a2v_(a2vec),
0015 iEtaLimits_(iEtaLimits),
0016 cut_(cut),
0017 minCharge_(minCharge),
0018 tFirst_(firstTimeSlice),
0019 tLast_(lastTimeSlice) {
0020 if (!validate())
0021 throw cms::Exception("Invalid HBHENegativeEFilter constructor arguments");
0022 }
0023
0024 bool HBHENegativeEFilter::validate() const {
0025 if (cut_.empty())
0026 return true;
0027
0028 const std::size_t nLimits(iEtaLimits_.size());
0029 if (nLimits >= static_cast<std::size_t>(UINT_MAX - 1U))
0030 return false;
0031 for (std::size_t i = 1; i < nLimits; ++i)
0032 if (!(iEtaLimits_[i - 1] < iEtaLimits_[i]))
0033 return false;
0034
0035 if (a1v_.size() != nLimits + 1)
0036 return false;
0037 if (a2v_.size() != nLimits + 1)
0038 return false;
0039
0040 const std::size_t sz = cut_.size();
0041 if (sz >= static_cast<std::size_t>(UINT_MAX - 1U))
0042 return false;
0043 for (std::size_t i = 1; i < sz; ++i)
0044 if (!(cut_[i - 1U].first < cut_[i].first))
0045 return false;
0046
0047 if (tFirst_ < 2U)
0048 return false;
0049 if (!(tFirst_ <= tLast_))
0050 return false;
0051
0052 return true;
0053 }
0054
0055 bool HBHENegativeEFilter::operator==(const HBHENegativeEFilter& r) const {
0056 if (cut_.empty() && r.cut_.empty())
0057 return true;
0058 else
0059 return a1v_ == r.a1v_ && a2v_ == r.a2v_ && iEtaLimits_ == r.iEtaLimits_ && cut_ == r.cut_ &&
0060 minCharge_ == r.minCharge_ && tFirst_ == r.tFirst_ && tLast_ == r.tLast_;
0061 }
0062
0063 unsigned HBHENegativeEFilter::getEtaIndex(const HcalDetId& id) const {
0064 const unsigned nLimits = iEtaLimits_.size();
0065 unsigned which(0U);
0066 if (nLimits) {
0067 const uint32_t uEta = std::abs(id.ieta());
0068 const uint32_t* limits(&iEtaLimits_[0]);
0069 for (; which < nLimits; ++which)
0070 if (uEta < limits[which])
0071 break;
0072 }
0073 return which;
0074 }
0075
0076 bool HBHENegativeEFilter::checkPassFilter(const HcalDetId& id, const double* ts, const unsigned lenTS) const {
0077 bool passes = true;
0078 const unsigned sz = cut_.size();
0079 if (sz) {
0080 double chargeInWindow = 0.0;
0081 for (unsigned i = tFirst_; i <= tLast_ && i < lenTS; ++i)
0082 chargeInWindow += ts[i];
0083 if (chargeInWindow >= minCharge_) {
0084
0085 const std::pair<double, double>* cut = &cut_[0];
0086 double cutValue = cut[0].second;
0087 if (sz > 1U) {
0088
0089 unsigned largerPoint = 0;
0090 for (; cut[largerPoint].first <= chargeInWindow; ++largerPoint) {
0091 }
0092
0093
0094 if (largerPoint >= sz)
0095 cutValue = cut[sz - 1U].second;
0096 else if (largerPoint) {
0097 const double slope = (cut[largerPoint].second - cut[largerPoint - 1U].second) /
0098 (cut[largerPoint].first - cut[largerPoint - 1U].first);
0099 cutValue = cut[largerPoint - 1U].second + slope * (chargeInWindow - cut[largerPoint - 1U].first);
0100 }
0101 }
0102
0103
0104 const unsigned itaIdx = getEtaIndex(id);
0105 const PiecewiseScalingPolynomial& a1(a1v_[itaIdx]);
0106 const PiecewiseScalingPolynomial& a2(a2v_[itaIdx]);
0107
0108 for (unsigned i = tFirst_; i <= tLast_ && i < lenTS && passes; ++i) {
0109 const double ecorr = ts[i] - a1(ts[i - 1U]) - a2(ts[i - 2U]);
0110 passes = ecorr >= cutValue;
0111 }
0112 }
0113 }
0114 return passes;
0115 }