HFSimpleTimeCheck.cc

Go to the documentation of this file.

#include <cstring>
#include <climits>

#include "DataFormats/HcalRecHit/interface/HcalSpecialTimes.h"

#include "RecoLocalCalo/HcalRecAlgos/interface/HFSimpleTimeCheck.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HFRecHitAuxSetter.h"

// Phase 1 rechit status bit assignments
#include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"

namespace {
  inline float build_rechit_time(const float weightedEnergySum,
                                 const float weightedSum,
                                 const float sum,
                                 const unsigned count,
                                 const float valueIfNothingWorks,
                                 bool* resultComesFromTDC) {
    if (weightedEnergySum > 0.f) {
      *resultComesFromTDC = true;
      return weightedSum / weightedEnergySum;
    } else if (count) {
      *resultComesFromTDC = true;
      return sum / count;
    } else {
      *resultComesFromTDC = false;
      return valueIfNothingWorks;
    }
  }
}  // namespace

HFSimpleTimeCheck::HFSimpleTimeCheck(const std::pair<float, float> tlimits[2],
                                     const float energyWeights[2 * HFAnodeStatus::N_POSSIBLE_STATES - 1][2],
                                     const unsigned i_soiPhase,
                                     const float i_timeShift,
                                     const float i_triseIfNoTDC,
                                     const float i_tfallIfNoTDC,
                                     const float i_minChargeForUndershoot,
                                     const float i_minChargeForOvershoot,
                                     const bool rejectAllFailures,
                                     const bool alwaysCalculateQAsymmetry)
    : soiPhase_(i_soiPhase),
      timeShift_(i_timeShift),
      triseIfNoTDC_(i_triseIfNoTDC),
      tfallIfNoTDC_(i_tfallIfNoTDC),
      minChargeForUndershoot_(i_minChargeForUndershoot),
      minChargeForOvershoot_(i_minChargeForOvershoot),
      rejectAllFailures_(rejectAllFailures),
      alwaysQAsym_(alwaysCalculateQAsymmetry) {
  tlimits_[0] = tlimits[0];
  tlimits_[1] = tlimits[1];
  float* to = &energyWeights_[0][0];
  const float* from = &energyWeights[0][0];
  memcpy(to, from, sizeof(energyWeights_));
}

unsigned HFSimpleTimeCheck::determineAnodeStatus(const unsigned ianode, const HFQIE10Info& anode, bool*) const {
  // Check if this anode has a dataframe error
  if (!anode.isDataframeOK())
    return HFAnodeStatus::HARDWARE_ERROR;

  // Check the time limits
  float trise = anode.timeRising();
  const bool timeIsKnown = !HcalSpecialTimes::isSpecial(trise);
  trise += timeShift_;
  if (timeIsKnown && tlimits_[ianode].first <= trise && trise <= tlimits_[ianode].second)
    return HFAnodeStatus::OK;
  else
    return HFAnodeStatus::FAILED_TIMING;
}

unsigned HFSimpleTimeCheck::mapStatusIntoIndex(const unsigned states[2]) const {
  unsigned eStates[2];
  eStates[0] = states[0];
  eStates[1] = states[1];
  if (!rejectAllFailures_)
    for (unsigned i = 0; i < 2; ++i)
      if (eStates[i] == HFAnodeStatus::FAILED_TIMING || eStates[i] == HFAnodeStatus::FAILED_OTHER)
        eStates[i] = HFAnodeStatus::OK;
  if (eStates[0] == HFAnodeStatus::OK)
    return eStates[1];
  else if (eStates[1] == HFAnodeStatus::OK)
    return HFAnodeStatus::N_POSSIBLE_STATES + eStates[0] - 1;
  else
    return UINT_MAX;
}

HFRecHit HFSimpleTimeCheck::reconstruct(const HFPreRecHit& prehit,
                                        const HcalCalibrations& /* calibs */,
                                        const bool flaggedBadInDB[2],
                                        const bool expectSingleAnodePMT) {
  HFRecHit rh;

  // Determine the status of each anode
  unsigned states[2] = {HFAnodeStatus::NOT_READ_OUT, HFAnodeStatus::NOT_READ_OUT};
  if (expectSingleAnodePMT)
    states[1] = HFAnodeStatus::NOT_DUAL;

  bool isTimingReliable[2] = {true, true};
  for (unsigned ianode = 0; ianode < 2; ++ianode) {
    if (flaggedBadInDB[ianode])
      states[ianode] = HFAnodeStatus::FLAGGED_BAD;
    else {
      const HFQIE10Info* anodeInfo = prehit.getHFQIE10Info(ianode);
      if (anodeInfo)
        states[ianode] = determineAnodeStatus(ianode, *anodeInfo, &isTimingReliable[ianode]);
    }
  }

  // Reconstruct energy and time
  const unsigned lookupInd = mapStatusIntoIndex(states);
  if (lookupInd != UINT_MAX) {
    // In this scope, at least one of states[i] is HFAnodeStatus::OK
    // or was mapped into that status by "mapStatusIntoIndex" method
    //
    const float* weights = &energyWeights_[lookupInd][0];
    float energy = 0.f, tfallWeightedEnergySum = 0.f, triseWeightedEnergySum = 0.f;
    float tfallWeightedSum = 0.f, triseWeightedSum = 0.f;
    float tfallSum = 0.f, triseSum = 0.f;
    unsigned tfallCount = 0, triseCount = 0;

    for (unsigned ianode = 0; ianode < 2; ++ianode) {
      const HFQIE10Info* anodeInfo = prehit.getHFQIE10Info(ianode);
      if (anodeInfo && weights[ianode] > 0.f) {
        const float weightedEnergy = weights[ianode] * anodeInfo->energy();
        energy += weightedEnergy;

        if (isTimingReliable[ianode] && states[ianode] != HFAnodeStatus::FAILED_TIMING) {
          float trise = anodeInfo->timeRising();
          if (!HcalSpecialTimes::isSpecial(trise)) {
            trise += timeShift_;
            triseSum += trise;
            ++triseCount;
            if (weightedEnergy > 0.f) {
              triseWeightedSum += trise * weightedEnergy;
              triseWeightedEnergySum += weightedEnergy;
            }
          }

          float tfall = anodeInfo->timeFalling();
          if (!HcalSpecialTimes::isSpecial(tfall)) {
            tfall += timeShift_;
            tfallSum += tfall;
            ++tfallCount;
            if (weightedEnergy > 0.f) {
              tfallWeightedSum += tfall * weightedEnergy;
              tfallWeightedEnergySum += weightedEnergy;
            }
          }
        }
      }
    }

    bool triseFromTDC = false;
    const float trise =
        build_rechit_time(triseWeightedEnergySum, triseWeightedSum, triseSum, triseCount, triseIfNoTDC_, &triseFromTDC);

    bool tfallFromTDC = false;
    const float tfall =
        build_rechit_time(tfallWeightedEnergySum, tfallWeightedSum, tfallSum, tfallCount, tfallIfNoTDC_, &tfallFromTDC);

    rh = HFRecHit(prehit.id(), energy, trise, tfall);
    HFRecHitAuxSetter::setAux(prehit, states, soiPhase_, &rh);

    // Set the "timing from TDC" flag
    const uint32_t flag = triseFromTDC ? 1U : 0U;
    rh.setFlagField(flag, HcalPhase1FlagLabels::TimingFromTDC);
  }

  return rh;
}