HFFlexibleTimeCheck.cc

Go to the documentation of this file.

#include <cmath>

#include "FWCore/Utilities/interface/Exception.h"
#include "DataFormats/HcalRecHit/interface/HcalSpecialTimes.h"
#include "DataFormats/HcalRecHit/interface/CaloRecHitAuxSetter.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HFRecHitAuxSetter.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HFFlexibleTimeCheck.h"

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

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

    // Require minimum charge for reliable timing measurement
    const float charge = anode.charge();
    const float minCharge = ianode ? pmtInfo_->minCharge1() :
                                     pmtInfo_->minCharge0();
    if (charge < minCharge)
    {
        *isTimingReliable = false;
        return HFAnodeStatus::OK;
    }

    // Special handling of under/overshoot TDC values, as well
    // as of DLL failures
    float trise = anode.timeRising();
    if ((trise == HcalSpecialTimes::UNKNOWN_T_UNDERSHOOT &&
         charge < minChargeForUndershoot()) ||
        (trise == HcalSpecialTimes::UNKNOWN_T_OVERSHOOT &&
         charge < minChargeForOvershoot()) ||
        trise == HcalSpecialTimes::UNKNOWN_T_DLL_FAILURE)
    {
        *isTimingReliable = false;
        return HFAnodeStatus::OK;
    }

    // Check if the rise time information is meaningful
    if (HcalSpecialTimes::isSpecial(trise))
        return HFAnodeStatus::FAILED_TIMING;

    // Figure out the timing cuts for this PMT
    const AbsHcalFunctor& minTimeShape = pmtInfo_->cut(
        ianode ? HFPhase1PMTData::T_1_MIN : HFPhase1PMTData::T_0_MIN);
    const AbsHcalFunctor& maxTimeShape = pmtInfo_->cut(
        ianode ? HFPhase1PMTData::T_1_MAX : HFPhase1PMTData::T_0_MAX);

    // Apply the timing cuts
    trise += timeShift();
    if (minTimeShape(charge) <= trise && trise <= maxTimeShape(charge))
        return HFAnodeStatus::OK;
    else
        return HFAnodeStatus::FAILED_TIMING;
}

HFRecHit HFFlexibleTimeCheck::reconstruct(const HFPreRecHit& prehit,
                                          const HcalCalibrations& calibs,
                                          const bool flaggedBadInDB[2],
                                          const bool expectSingleAnodePMT)
{
    // The algorithm must be configured by now
    if (!algoConf_)
        throw cms::Exception("HFPhase1BadConfig")
            << "In HFFlexibleTimeCheck::reconstruct: algorithm is not configured";

    // Fetch the algorithm configuration data for this PMT
    pmtInfo_ = &algoConf_->at(prehit.id());

    // Run the reconstruction algorithm from the base class
    HFRecHit rh = HFSimpleTimeCheck::reconstruct(
        prehit, calibs, flaggedBadInDB, expectSingleAnodePMT);

    if (rh.id().rawId())
    {
        // Check the charge asymmetry between the two anodes
        bool setAsymmetryFlag = true;
        if (!alwaysCalculatingQAsym())
        {
            using namespace CaloRecHitAuxSetter;

            const unsigned st0 = getField(rh.aux(), HFRecHitAuxSetter::MASK_STATUS,
                                          HFRecHitAuxSetter::OFF_STATUS);
            const unsigned st1 = getField(rh.getAuxHF(), HFRecHitAuxSetter::MASK_STATUS,
                                          HFRecHitAuxSetter::OFF_STATUS);
            setAsymmetryFlag = st0 == HFAnodeStatus::OK && st1 == HFAnodeStatus::OK;
        }

        if (setAsymmetryFlag)
        {
            bool passesAsymmetryCut = true;
            const std::pair<float,bool> qAsymm =
                prehit.chargeAsymmetry(pmtInfo_->minChargeAsymm());
            if (qAsymm.second)
            {
                const float q = prehit.charge();
                const float minAsymm = (pmtInfo_->cut(HFPhase1PMTData::ASYMM_MIN))(q);
                const float maxAsymm = (pmtInfo_->cut(HFPhase1PMTData::ASYMM_MAX))(q);
                passesAsymmetryCut = minAsymm <= qAsymm.first && qAsymm.first <= maxAsymm;
            }
            if (!passesAsymmetryCut)
                rh.setFlagField(1U, HcalPhase1FlagLabels::HFSignalAsymmetry);
        }
    }

    return rh;
}