|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
#include <cmath>
#include <climits>
#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h"
HBHENegativeEFilter::HBHENegativeEFilter(const std::vector<PiecewiseScalingPolynomial>& a1vec,
const std::vector<PiecewiseScalingPolynomial>& a2vec,
const std::vector<uint32_t>& iEtaLimits,
const std::vector<std::pair<double, double> >& cut,
const double minCharge,
const unsigned firstTimeSlice,
const unsigned lastTimeSlice)
: a1v_(a1vec),
a2v_(a2vec),
iEtaLimits_(iEtaLimits),
cut_(cut),
minCharge_(minCharge),
tFirst_(firstTimeSlice),
tLast_(lastTimeSlice) {
if (!validate())
throw cms::Exception("Invalid HBHENegativeEFilter constructor arguments");
}
bool HBHENegativeEFilter::validate() const {
if (cut_.empty())
return true;
const std::size_t nLimits(iEtaLimits_.size());
if (nLimits >= static_cast<std::size_t>(UINT_MAX - 1U))
return false;
for (std::size_t i = 1; i < nLimits; ++i)
if (!(iEtaLimits_[i - 1] < iEtaLimits_[i]))
return false;
if (a1v_.size() != nLimits + 1)
return false;
if (a2v_.size() != nLimits + 1)
return false;
const std::size_t sz = cut_.size();
if (sz >= static_cast<std::size_t>(UINT_MAX - 1U))
return false;
for (std::size_t i = 1; i < sz; ++i)
if (!(cut_[i - 1U].first < cut_[i].first))
return false;
if (tFirst_ < 2U)
return false;
if (!(tFirst_ <= tLast_))
return false;
return true;
}
bool HBHENegativeEFilter::operator==(const HBHENegativeEFilter& r) const {
if (cut_.empty() && r.cut_.empty())
return true;
else
return a1v_ == r.a1v_ && a2v_ == r.a2v_ && iEtaLimits_ == r.iEtaLimits_ && cut_ == r.cut_ &&
minCharge_ == r.minCharge_ && tFirst_ == r.tFirst_ && tLast_ == r.tLast_;
}
unsigned HBHENegativeEFilter::getEtaIndex(const HcalDetId& id) const {
const unsigned nLimits = iEtaLimits_.size();
unsigned which(0U);
if (nLimits) {
const uint32_t uEta = std::abs(id.ieta());
const uint32_t* limits(&iEtaLimits_[0]);
for (; which < nLimits; ++which)
if (uEta < limits[which])
break;
}
return which;
}
bool HBHENegativeEFilter::checkPassFilter(const HcalDetId& id, const double* ts, const unsigned lenTS) const {
bool passes = true;
const unsigned sz = cut_.size();
if (sz) {
double chargeInWindow = 0.0;
for (unsigned i = tFirst_; i <= tLast_ && i < lenTS; ++i)
chargeInWindow += ts[i];
if (chargeInWindow >= minCharge_) {
// Figure out the cut value for this charge
const std::pair<double, double>* cut = &cut_[0];
double cutValue = cut[0].second;
if (sz > 1U) {
// First point larger than charge
unsigned largerPoint = 0;
for (; cut[largerPoint].first <= chargeInWindow; ++largerPoint) {
}
// Constant extrapolation beyond min and max coords
if (largerPoint >= sz)
cutValue = cut[sz - 1U].second;
else if (largerPoint) {
const double slope = (cut[largerPoint].second - cut[largerPoint - 1U].second) /
(cut[largerPoint].first - cut[largerPoint - 1U].first);
cutValue = cut[largerPoint - 1U].second + slope * (chargeInWindow - cut[largerPoint - 1U].first);
}
}
// Compare the modified time slices with the cut
const unsigned itaIdx = getEtaIndex(id);
const PiecewiseScalingPolynomial& a1(a1v_[itaIdx]);
const PiecewiseScalingPolynomial& a2(a2v_[itaIdx]);
for (unsigned i = tFirst_; i <= tLast_ && i < lenTS && passes; ++i) {
const double ecorr = ts[i] - a1(ts[i - 1U]) - a2(ts[i - 2U]);
passes = ecorr >= cutValue;
}
}
}
return passes;
}
|