HcalRecAlgos/src/HFFlexibleTimeCheck.cc

0001 #include <cmath>
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "DataFormats/HcalRecHit/interface/HcalSpecialTimes.h"
0005 #include "DataFormats/HcalRecHit/interface/CaloRecHitAuxSetter.h"
0006 #include "RecoLocalCalo/HcalRecAlgos/interface/HFRecHitAuxSetter.h"
0007 #include "RecoLocalCalo/HcalRecAlgos/interface/HFFlexibleTimeCheck.h"
0008
0009 // Phase 1 rechit status bit assignments
0010 #include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"
0011
0012 unsigned HFFlexibleTimeCheck::determineAnodeStatus(const unsigned ianode,
0013                                                    const HFQIE10Info& anode,
0014                                                    bool* isTimingReliable) const {
0015   // Return quickly if this anode has a dataframe error
0016   if (!anode.isDataframeOK())
0017     return HFAnodeStatus::HARDWARE_ERROR;
0018
0019   // Require minimum charge for reliable timing measurement
0020   const float charge = anode.charge();
0021   const float minCharge = ianode ? pmtInfo_->minCharge1() : pmtInfo_->minCharge0();
0022   if (charge < minCharge) {
0023     *isTimingReliable = false;
0024     return HFAnodeStatus::OK;
0025   }
0026
0027   // Special handling of under/overshoot TDC values, as well
0028   // as of DLL failures
0029   float trise = anode.timeRising();
0030   if ((trise == HcalSpecialTimes::UNKNOWN_T_UNDERSHOOT && charge < minChargeForUndershoot()) ||
0031       (trise == HcalSpecialTimes::UNKNOWN_T_OVERSHOOT && charge < minChargeForOvershoot()) ||
0032       trise == HcalSpecialTimes::UNKNOWN_T_DLL_FAILURE) {
0033     *isTimingReliable = false;
0034     return HFAnodeStatus::OK;
0035   }
0036
0037   // Check if the rise time information is meaningful
0038   if (HcalSpecialTimes::isSpecial(trise))
0039     return HFAnodeStatus::FAILED_TIMING;
0040
0041   // Figure out the timing cuts for this PMT
0042   const AbsHcalFunctor& minTimeShape = pmtInfo_->cut(ianode ? HFPhase1PMTData::T_1_MIN : HFPhase1PMTData::T_0_MIN);
0043   const AbsHcalFunctor& maxTimeShape = pmtInfo_->cut(ianode ? HFPhase1PMTData::T_1_MAX : HFPhase1PMTData::T_0_MAX);
0044
0045   // Apply the timing cuts
0046   trise += timeShift();
0047   if (minTimeShape(charge) <= trise && trise <= maxTimeShape(charge))
0048     return HFAnodeStatus::OK;
0049   else
0050     return HFAnodeStatus::FAILED_TIMING;
0051 }
0052
0053 HFRecHit HFFlexibleTimeCheck::reconstruct(const HFPreRecHit& prehit,
0054                                           const HcalCalibrations& calibs,
0055                                           const bool flaggedBadInDB[2],
0056                                           const bool expectSingleAnodePMT) {
0057   // The algorithm must be configured by now
0058   if (!algoConf_)
0059     throw cms::Exception("HFPhase1BadConfig") << "In HFFlexibleTimeCheck::reconstruct: algorithm is not configured";
0060
0061   // Fetch the algorithm configuration data for this PMT
0062   pmtInfo_ = &algoConf_->at(prehit.id());
0063
0064   // Run the reconstruction algorithm from the base class
0065   HFRecHit rh = HFSimpleTimeCheck::reconstruct(prehit, calibs, flaggedBadInDB, expectSingleAnodePMT);
0066
0067   if (rh.id().rawId()) {
0068     // Check the charge asymmetry between the two anodes
0069     bool setAsymmetryFlag = true;
0070     if (!alwaysCalculatingQAsym()) {
0071       using namespace CaloRecHitAuxSetter;
0072
0073       const unsigned st0 = getField(rh.aux(), HFRecHitAuxSetter::MASK_STATUS, HFRecHitAuxSetter::OFF_STATUS);
0074       const unsigned st1 = getField(rh.getAuxHF(), HFRecHitAuxSetter::MASK_STATUS, HFRecHitAuxSetter::OFF_STATUS);
0075       setAsymmetryFlag = st0 == HFAnodeStatus::OK && st1 == HFAnodeStatus::OK;
0076     }
0077
0078     if (setAsymmetryFlag) {
0079       bool passesAsymmetryCut = true;
0080       const std::pair<float, bool> qAsymm = prehit.chargeAsymmetry(pmtInfo_->minChargeAsymm());
0081       if (qAsymm.second) {
0082         const float q = prehit.charge();
0083         const float minAsymm = (pmtInfo_->cut(HFPhase1PMTData::ASYMM_MIN))(q);
0084         const float maxAsymm = (pmtInfo_->cut(HFPhase1PMTData::ASYMM_MAX))(q);
0085         passesAsymmetryCut = minAsymm <= qAsymm.first && qAsymm.first <= maxAsymm;
0086       }
0087       if (!passesAsymmetryCut)
0088         rh.setFlagField(1U, HcalPhase1FlagLabels::HFSignalAsymmetry);
0089     }
0090   }
0091
0092   return rh;
0093 }