File indexing completed on 2023-03-17 11:18:52
0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HFTimingTrustFlag.h"
0002 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0003
0004 namespace HFTimingTrust {
0005
0006
0007
0008 double timerr_hf(double rectime, double ampl);
0009
0010 template <class T, class V>
0011 void checkHFTimErr(T& rechit, const V& digi, int level1, int level2) {
0012
0013 double rectime = rechit.time();
0014
0015 if (rectime > -100 && rectime < 250) {
0016
0017 double ampl = 0;
0018 double sum = 0;
0019 int maxI = -1;
0020 for (int i = 0; i < digi.size(); ++i) {
0021 sum += digi.sample(i).nominal_fC();
0022 if (digi.sample(i).nominal_fC() > ampl) {
0023 ampl = digi.sample(i).nominal_fC();
0024 maxI = i;
0025 }
0026 }
0027 if (ampl > 1 && maxI > 0 && maxI < digi.size() - 1) {
0028 ampl = ampl + digi.sample(maxI - 1).nominal_fC() + digi.sample(maxI + 1).nominal_fC();
0029 ampl -= (sum - ampl) * 3.0 / (digi.size() - 3);
0030 if (ampl > 3) {
0031 int timerr = (int)timerr_hf(rectime, ampl);
0032 uint32_t status = 0;
0033 if (timerr < 0)
0034 status = 3;
0035 else if (timerr <= level1)
0036 status = 0;
0037 else if (timerr <= level2)
0038 status = 1;
0039 else
0040 status = 2;
0041 rechit.setFlagField(status, HcalCaloFlagLabels::HFTimingTrustBits, 2);
0042 return;
0043 }
0044 }
0045
0046 }
0047
0048 rechit.setFlagField(3, HcalCaloFlagLabels::HFTimingTrustBits, 2);
0049 return;
0050 }
0051
0052 static const float hfterrlut[195] = {
0053 3.42, 3.04, 9.18, 8.97, 8.49, 8.08, 8.01, 8.30, 8.75, 8.22, 7.21, 5.04, 2.98, 2.04, 1.22, 7.58, 7.79, 7.11,
0054 6.93, 7.03, 7.27, 7.23, 6.53, 4.59, 2.40, 1.46, 1.31, 0.42, 5.95, 6.48, 6.29, 5.84, 5.97, 6.31, 6.00, 4.37,
0055 2.37, 1.03, 0.72, 0.81, 0.27, 3.98, 5.57, 5.04, 5.10, 5.21, 5.18, 4.22, 2.23, 1.07, 0.66, 0.40, 0.48, 0.17,
0056 2.51, 4.70, 4.28, 4.29, 4.36, 3.84, 2.40, 1.15, 0.68, 0.40, 0.24, 0.29, 0.11, 0.81, 3.71, 3.47, 3.48, 3.52,
0057 2.58, 1.25, 0.71, 0.41, 0.26, 0.16, 0.16, 0.08, 0.27, 2.88, 2.63, 2.76, 2.33, 1.31, 0.72, 0.44, 0.27, 0.16,
0058 0.11, 0.10, 0.06, 0.15, 2.11, 2.00, 1.84, 1.46, 0.79, 0.45, 0.26, 0.17, 0.10, 0.08, 0.05, 0.04, 0.10, 1.58,
0059 1.49, 1.25, 0.90, 0.48, 0.29, 0.17, 0.10, 0.06, 0.06, 0.02, 0.03, 0.06, 1.26, 1.03, 0.77, 0.57, 0.30, 0.18,
0060 0.11, 0.06, 0.04, 0.05, 0.01, 0.02, 0.04, 0.98, 0.66, 0.47, 0.39, 0.18, 0.11, 0.07, 0.04, 0.03, 0.04, 0.01,
0061 0.02, 0.02, 0.86, 0.44, 0.30, 0.27, 0.11, 0.07, 0.04, 0.03, 0.02, 0.04, 0.01, 0.02, 0.02, 0.80, 0.30, 0.21,
0062 0.17, 0.07, 0.04, 0.03, 0.02, 0.01, 0.04, 0.01, 0.02, 0.01, 0.76, 0.22, 0.17, 0.12, 0.05, 0.03, 0.02, 0.01,
0063 0.01, 0.04, 0.01, 0.02, 0.01, 0.76, 0.17, 0.14, 0.09, 0.03, 0.02, 0.01, 0.01, 0.01, 0.04};
0064
0065 double timerr_hf(double rectime, double ampl) {
0066 int itim, iampl, index;
0067 double tim;
0068 if (rectime > 0)
0069 tim = rectime - ((int)(rectime / 25)) * 25;
0070 else
0071 tim = rectime - ((int)(rectime / 25)) * 25 + 25;
0072 itim = (int)(tim / 2.0);
0073
0074 iampl = 0;
0075 static const double bampl[15] = {3, 5, 8, 12, 19, 30, 47, 73, 115, 182, 287, 452, 712, 1120, 1766};
0076 if (ampl >= bampl[14])
0077 iampl = 14;
0078 else {
0079 for (int i = 1; i <= 14; i++) {
0080 if (ampl < bampl[i]) {
0081 iampl = i - 1;
0082 break;
0083 }
0084 }
0085 }
0086
0087 index = itim + iampl * 13;
0088
0089 double y1 = hfterrlut[index];
0090 double y2 = 0;
0091 double v1 = y1;
0092 if (itim < 12) {
0093 y2 = hfterrlut[index + 1];
0094 v1 = y1 + (y2 - y1) * (tim / 2.0 - (float)itim);
0095 }
0096 double yy1 = 0;
0097 double yy2 = 0;
0098 double v2 = 0;
0099 if (iampl < 14) {
0100 yy1 = hfterrlut[index + 13];
0101 if (itim == 12)
0102 v2 = yy1;
0103 else {
0104 yy2 = hfterrlut[index + 14];
0105 v2 = yy1 + (yy2 - yy1) * (tim / 2.0 - (float)itim);
0106 }
0107 v1 = v1 + (v2 - v1) * (ampl - bampl[iampl]) / (bampl[iampl + 1] - bampl[iampl]);
0108 }
0109
0110 return v1;
0111 }
0112
0113 }
0114
0115 using namespace HFTimingTrust;
0116
0117 HFTimingTrustFlag::HFTimingTrustFlag() {
0118 HFTimingTrustLevel1_ = 1;
0119 HFTimingTrustLevel2_ = 4;
0120 }
0121
0122 HFTimingTrustFlag::HFTimingTrustFlag(int level1, int level2) {
0123 HFTimingTrustLevel1_ = level1;
0124 HFTimingTrustLevel2_ = level2;
0125 }
0126
0127 HFTimingTrustFlag::~HFTimingTrustFlag() {}
0128
0129 void HFTimingTrustFlag::setHFTimingTrustFlag(HFRecHit& rechit, const HFDataFrame& digi) {
0130 checkHFTimErr<HFRecHit, HFDataFrame>(rechit, digi, HFTimingTrustLevel1_, HFTimingTrustLevel2_);
0131 return;
0132 }