HcalObjects/src/HBHENegativeEFilter.cc

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       // Figure out the cut value for this charge
0085       const std::pair<double, double>* cut = &cut_[0];
0086       double cutValue = cut[0].second;
0087       if (sz > 1U) {
0088         // First point larger than charge
0089         unsigned largerPoint = 0;
0090         for (; cut[largerPoint].first <= chargeInWindow; ++largerPoint) {
0091         }
0092
0093         // Constant extrapolation beyond min and max coords
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       // Compare the modified time slices with the cut
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 }