HcalRecAlgos/src/HFSimpleTimeCheck.cc

0001 #include <cstring>
0002 #include <climits>
0003
0004 #include "DataFormats/HcalRecHit/interface/HcalSpecialTimes.h"
0005
0006 #include "RecoLocalCalo/HcalRecAlgos/interface/HFSimpleTimeCheck.h"
0007 #include "RecoLocalCalo/HcalRecAlgos/interface/HFRecHitAuxSetter.h"
0008
0009 // Phase 1 rechit status bit assignments
0010 #include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"
0011
0012 namespace {
0013   inline float build_rechit_time(const float weightedEnergySum,
0014                                  const float weightedSum,
0015                                  const float sum,
0016                                  const unsigned count,
0017                                  const float valueIfNothingWorks,
0018                                  bool* resultComesFromTDC) {
0019     if (weightedEnergySum > 0.f) {
0020       *resultComesFromTDC = true;
0021       return weightedSum / weightedEnergySum;
0022     } else if (count) {
0023       *resultComesFromTDC = true;
0024       return sum / count;
0025     } else {
0026       *resultComesFromTDC = false;
0027       return valueIfNothingWorks;
0028     }
0029   }
0030 }  // namespace
0031
0032 HFSimpleTimeCheck::HFSimpleTimeCheck(const std::pair<float, float> tlimits[2],
0033                                      const float energyWeights[2 * HFAnodeStatus::N_POSSIBLE_STATES - 1][2],
0034                                      const unsigned i_soiPhase,
0035                                      const float i_timeShift,
0036                                      const float i_triseIfNoTDC,
0037                                      const float i_tfallIfNoTDC,
0038                                      const float i_minChargeForUndershoot,
0039                                      const float i_minChargeForOvershoot,
0040                                      const bool rejectAllFailures,
0041                                      const bool alwaysCalculateQAsymmetry)
0042     : soiPhase_(i_soiPhase),
0043       timeShift_(i_timeShift),
0044       triseIfNoTDC_(i_triseIfNoTDC),
0045       tfallIfNoTDC_(i_tfallIfNoTDC),
0046       minChargeForUndershoot_(i_minChargeForUndershoot),
0047       minChargeForOvershoot_(i_minChargeForOvershoot),
0048       rejectAllFailures_(rejectAllFailures),
0049       alwaysQAsym_(alwaysCalculateQAsymmetry) {
0050   tlimits_[0] = tlimits[0];
0051   tlimits_[1] = tlimits[1];
0052   float* to = &energyWeights_[0][0];
0053   const float* from = &energyWeights[0][0];
0054   memcpy(to, from, sizeof(energyWeights_));
0055 }
0056
0057 unsigned HFSimpleTimeCheck::determineAnodeStatus(const unsigned ianode, const HFQIE10Info& anode, bool*) const {
0058   // Check if this anode has a dataframe error
0059   if (!anode.isDataframeOK())
0060     return HFAnodeStatus::HARDWARE_ERROR;
0061
0062   // Check the time limits
0063   float trise = anode.timeRising();
0064   const bool timeIsKnown = !HcalSpecialTimes::isSpecial(trise);
0065   trise += timeShift_;
0066   if (timeIsKnown && tlimits_[ianode].first <= trise && trise <= tlimits_[ianode].second)
0067     return HFAnodeStatus::OK;
0068   else
0069     return HFAnodeStatus::FAILED_TIMING;
0070 }
0071
0072 unsigned HFSimpleTimeCheck::mapStatusIntoIndex(const unsigned states[2]) const {
0073   unsigned eStates[2];
0074   eStates[0] = states[0];
0075   eStates[1] = states[1];
0076   if (!rejectAllFailures_)
0077     for (unsigned i = 0; i < 2; ++i)
0078       if (eStates[i] == HFAnodeStatus::FAILED_TIMING || eStates[i] == HFAnodeStatus::FAILED_OTHER)
0079         eStates[i] = HFAnodeStatus::OK;
0080   if (eStates[0] == HFAnodeStatus::OK)
0081     return eStates[1];
0082   else if (eStates[1] == HFAnodeStatus::OK)
0083     return HFAnodeStatus::N_POSSIBLE_STATES + eStates[0] - 1;
0084   else
0085     return UINT_MAX;
0086 }
0087
0088 HFRecHit HFSimpleTimeCheck::reconstruct(const HFPreRecHit& prehit,
0089                                         const HcalCalibrations& /* calibs */,
0090                                         const bool flaggedBadInDB[2],
0091                                         const bool expectSingleAnodePMT) {
0092   HFRecHit rh;
0093
0094   // Determine the status of each anode
0095   unsigned states[2] = {HFAnodeStatus::NOT_READ_OUT, HFAnodeStatus::NOT_READ_OUT};
0096   if (expectSingleAnodePMT)
0097     states[1] = HFAnodeStatus::NOT_DUAL;
0098
0099   bool isTimingReliable[2] = {true, true};
0100   for (unsigned ianode = 0; ianode < 2; ++ianode) {
0101     if (flaggedBadInDB[ianode])
0102       states[ianode] = HFAnodeStatus::FLAGGED_BAD;
0103     else {
0104       const HFQIE10Info* anodeInfo = prehit.getHFQIE10Info(ianode);
0105       if (anodeInfo)
0106         states[ianode] = determineAnodeStatus(ianode, *anodeInfo, &isTimingReliable[ianode]);
0107     }
0108   }
0109
0110   // Reconstruct energy and time
0111   const unsigned lookupInd = mapStatusIntoIndex(states);
0112   if (lookupInd != UINT_MAX) {
0113     // In this scope, at least one of states[i] is HFAnodeStatus::OK
0114     // or was mapped into that status by "mapStatusIntoIndex" method
0115     //
0116     const float* weights = &energyWeights_[lookupInd][0];
0117     float energy = 0.f, tfallWeightedEnergySum = 0.f, triseWeightedEnergySum = 0.f;
0118     float tfallWeightedSum = 0.f, triseWeightedSum = 0.f;
0119     float tfallSum = 0.f, triseSum = 0.f;
0120     unsigned tfallCount = 0, triseCount = 0;
0121
0122     for (unsigned ianode = 0; ianode < 2; ++ianode) {
0123       const HFQIE10Info* anodeInfo = prehit.getHFQIE10Info(ianode);
0124       if (anodeInfo && weights[ianode] > 0.f) {
0125         const float weightedEnergy = weights[ianode] * anodeInfo->energy();
0126         energy += weightedEnergy;
0127
0128         if (isTimingReliable[ianode] && states[ianode] != HFAnodeStatus::FAILED_TIMING) {
0129           float trise = anodeInfo->timeRising();
0130           if (!HcalSpecialTimes::isSpecial(trise)) {
0131             trise += timeShift_;
0132             triseSum += trise;
0133             ++triseCount;
0134             if (weightedEnergy > 0.f) {
0135               triseWeightedSum += trise * weightedEnergy;
0136               triseWeightedEnergySum += weightedEnergy;
0137             }
0138           }
0139
0140           float tfall = anodeInfo->timeFalling();
0141           if (!HcalSpecialTimes::isSpecial(tfall)) {
0142             tfall += timeShift_;
0143             tfallSum += tfall;
0144             ++tfallCount;
0145             if (weightedEnergy > 0.f) {
0146               tfallWeightedSum += tfall * weightedEnergy;
0147               tfallWeightedEnergySum += weightedEnergy;
0148             }
0149           }
0150         }
0151       }
0152     }
0153
0154     bool triseFromTDC = false;
0155     const float trise =
0156         build_rechit_time(triseWeightedEnergySum, triseWeightedSum, triseSum, triseCount, triseIfNoTDC_, &triseFromTDC);
0157
0158     bool tfallFromTDC = false;
0159     const float tfall =
0160         build_rechit_time(tfallWeightedEnergySum, tfallWeightedSum, tfallSum, tfallCount, tfallIfNoTDC_, &tfallFromTDC);
0161
0162     rh = HFRecHit(prehit.id(), energy, trise, tfall);
0163     HFRecHitAuxSetter::setAux(prehit, states, soiPhase_, &rh);
0164
0165     // Set the "timing from TDC" flag
0166     const uint32_t flag = triseFromTDC ? 1U : 0U;
0167     rh.setFlagField(flag, HcalPhase1FlagLabels::TimingFromTDC);
0168   }
0169
0170   return rh;
0171 }