HcalRecAlgos/src/HFTimingTrustFlag.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HFTimingTrustFlag.h"
0002 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0003
0004 namespace HFTimingTrust {
0005   // Template class checks HF timing, sets rechit
0006   // timing status bits according to its uncertainty;
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     // Get rechit timing
0013     double rectime = rechit.time();
0014
0015     if (rectime > -100 && rectime < 250) {
0016       // Get signal from digi
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;  // unreconstructable time
0035           else if (timerr <= level1)
0036             status = 0;  // time reconstructed; precision better than level1 value
0037           else if (timerr <= level2)
0038             status = 1;  // precision worse than level1 value
0039           else
0040             status = 2;  //precision worse than level 2 value
0041           rechit.setFlagField(status, HcalCaloFlagLabels::HFTimingTrustBits, 2);
0042           return;
0043         }
0044       }
0045
0046     }  // if (rectime > -100 && rectime < -250)
0047     // if rectime outside expected range, set flag field to 3 (unreconstructable)?
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 }  // namespace HFTimingTrust
0114
0115 using namespace HFTimingTrust;
0116
0117 HFTimingTrustFlag::HFTimingTrustFlag() {
0118   HFTimingTrustLevel1_ = 1;  // time precision 1ns
0119   HFTimingTrustLevel2_ = 4;  // time precision 4ns
0120 }
0121
0122 HFTimingTrustFlag::HFTimingTrustFlag(int level1, int level2) {
0123   HFTimingTrustLevel1_ = level1;  // allow user to set t-trust level
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 }