File indexing completed on 2023-03-17 11:18:53
0001 #include <cassert>
0002 #include <algorithm>
0003
0004 #include "DataFormats/HcalRecHit/interface/HcalSpecialTimes.h"
0005 #include "DataFormats/HcalRecHit/interface/HBHERecHitAuxSetter.h"
0006 #include "DataFormats/HcalRecHit/interface/CaloRecHitAuxSetter.h"
0007
0008 #include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"
0009
0010 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0011
0012 #include "RecoLocalCalo/HcalRecAlgos/interface/SimplePlan1RechitCombiner.h"
0013
0014 SimplePlan1RechitCombiner::SimplePlan1RechitCombiner() : topo_(nullptr) {}
0015
0016 void SimplePlan1RechitCombiner::setTopo(const HcalTopology* topo) { topo_ = topo; }
0017
0018 void SimplePlan1RechitCombiner::clear() { rechitMap_.clear(); }
0019
0020 void SimplePlan1RechitCombiner::add(const HBHERecHit& rh) {
0021 if (!CaloRecHitAuxSetter::getBit(rh.auxPhase1(), HBHERecHitAuxSetter::OFF_DROPPED))
0022 rechitMap_.push_back(MapItem(mapRechit(rh.id()), &rh));
0023 }
0024
0025 void SimplePlan1RechitCombiner::combine(HBHERecHitCollection* toFill) {
0026 if (!rechitMap_.empty()) {
0027 std::sort(rechitMap_.begin(), rechitMap_.end());
0028
0029 HcalDetId oldId(rechitMap_[0].first);
0030 ptrbuf_.clear();
0031 ptrbuf_.push_back(rechitMap_[0].second);
0032
0033 const std::size_t nInput = rechitMap_.size();
0034 for (std::size_t i = 1; i < nInput; ++i) {
0035 if (rechitMap_[i].first != oldId) {
0036 const HBHERecHit& rh = makeRechit(oldId, ptrbuf_);
0037 if (rh.id().rawId())
0038 toFill->push_back(rh);
0039 oldId = rechitMap_[i].first;
0040 ptrbuf_.clear();
0041 }
0042 ptrbuf_.push_back(rechitMap_[i].second);
0043 }
0044
0045 const HBHERecHit& rh = makeRechit(oldId, ptrbuf_);
0046 if (rh.id().rawId())
0047 toFill->push_back(rh);
0048 }
0049 }
0050
0051 HcalDetId SimplePlan1RechitCombiner::mapRechit(HcalDetId from) const { return topo_->mergedDepthDetId(from); }
0052
0053 HBHERecHit SimplePlan1RechitCombiner::makeRechit(const HcalDetId idToMake,
0054 const std::vector<const HBHERecHit*>& rechits) const {
0055 constexpr unsigned MAXLEN = 8U;
0056 constexpr float TIME_IF_NO_ENERGY = -999.f;
0057
0058 const unsigned nRecHits = rechits.size();
0059 assert(nRecHits);
0060 assert(nRecHits <= MAXLEN);
0061
0062
0063 double energy = 0.0, eraw = 0.0, eaux = 0.0, chisq = 0.0;
0064 FPair times[MAXLEN], adctimes[MAXLEN];
0065 unsigned nADCTimes = 0;
0066
0067 for (unsigned i = 0; i < nRecHits; ++i) {
0068 const HBHERecHit& rh(*rechits[i]);
0069 const float e = rh.energy();
0070 energy += e;
0071 eraw += rh.eraw();
0072 eaux += rh.eaux();
0073 chisq += rh.chi2();
0074 times[i].first = rh.time();
0075 times[i].second = e;
0076
0077 const float tADC = rh.timeFalling();
0078 if (!HcalSpecialTimes::isSpecial(tADC)) {
0079 adctimes[nADCTimes].first = tADC;
0080 adctimes[nADCTimes].second = e;
0081 ++nADCTimes;
0082 }
0083 }
0084
0085 HBHERecHit rh(idToMake,
0086 energy,
0087 energyWeightedAverage(times, nRecHits, TIME_IF_NO_ENERGY),
0088 energyWeightedAverage(adctimes, nADCTimes, HcalSpecialTimes::UNKNOWN_T_NOTDC));
0089 rh.setRawEnergy(eraw);
0090 rh.setAuxEnergy(eaux);
0091 rh.setChiSquared(chisq);
0092
0093
0094 combineAuxInfo(rechits, &rh);
0095
0096 return rh;
0097 }
0098
0099 void SimplePlan1RechitCombiner::combineAuxInfo(const std::vector<const HBHERecHit*>& rechits, HBHERecHit* rh) const {
0100 using namespace CaloRecHitAuxSetter;
0101 using namespace HcalPhase1FlagLabels;
0102
0103
0104
0105 constexpr unsigned MAXLEN = 8U;
0106
0107 const unsigned nRecHits = rechits.size();
0108 assert(nRecHits);
0109 assert(nRecHits <= MAXLEN);
0110
0111 uint32_t flags = 0, auxPhase1 = 0;
0112 unsigned tripleFitCount = 0;
0113 unsigned soiVote[HBHERecHitAuxSetter::MASK_SOI + 1U] = {0};
0114 unsigned capidVote[HBHERecHitAuxSetter::MASK_CAPID + 1U] = {0};
0115
0116
0117 for (unsigned i = 0; i < nRecHits; ++i) {
0118 const HBHERecHit& rh(*rechits[i]);
0119 const uint32_t rhflags = rh.flags();
0120 const uint32_t rhAuxPhase1 = rh.auxPhase1();
0121
0122 orBit(&auxPhase1, HBHERecHitAuxSetter::OFF_LINK_ERR, getBit(rhAuxPhase1, HBHERecHitAuxSetter::OFF_LINK_ERR));
0123 orBit(&auxPhase1, HBHERecHitAuxSetter::OFF_CAPID_ERR, getBit(rhAuxPhase1, HBHERecHitAuxSetter::OFF_CAPID_ERR));
0124
0125 const unsigned soi = getField(rhAuxPhase1, HBHERecHitAuxSetter::MASK_SOI, HBHERecHitAuxSetter::OFF_SOI);
0126 soiVote[soi]++;
0127
0128 const unsigned capid = getField(rhAuxPhase1, HBHERecHitAuxSetter::MASK_CAPID, HBHERecHitAuxSetter::OFF_CAPID);
0129 capidVote[capid]++;
0130
0131 if (getBit(rhflags, HBHEStatusFlag::HBHEPulseFitBit))
0132 ++tripleFitCount;
0133
0134
0135
0136 flags |= rhflags;
0137 }
0138
0139 unsigned* pmaxsoi = std::max_element(soiVote, soiVote + sizeof(soiVote) / sizeof(soiVote[0]));
0140 const unsigned soi = std::distance(&soiVote[0], pmaxsoi);
0141 setField(&auxPhase1, HBHERecHitAuxSetter::MASK_SOI, HBHERecHitAuxSetter::OFF_SOI, soi);
0142
0143 unsigned* pmaxcapid = std::max_element(capidVote, capidVote + sizeof(capidVote) / sizeof(capidVote[0]));
0144 const unsigned capid = std::distance(&capidVote[0], pmaxcapid);
0145 setField(&auxPhase1, HBHERecHitAuxSetter::MASK_CAPID, HBHERecHitAuxSetter::OFF_CAPID, capid);
0146
0147
0148 setField(&auxPhase1, HBHERecHitAuxSetter::MASK_NSAMPLES, HBHERecHitAuxSetter::OFF_ADC, tripleFitCount);
0149
0150
0151 setField(&auxPhase1, HBHERecHitAuxSetter::MASK_NSAMPLES, HBHERecHitAuxSetter::OFF_NSAMPLES, nRecHits);
0152
0153
0154 setBit(&auxPhase1, HBHERecHitAuxSetter::OFF_TDC_TIME, true);
0155
0156
0157 setBit(&auxPhase1, HBHERecHitAuxSetter::OFF_COMBINED, true);
0158
0159
0160 rh->setFlags(flags);
0161 rh->setAuxPhase1(auxPhase1);
0162
0163
0164
0165 unsigned depthValues[MAXLEN];
0166 for (unsigned i = 0; i < nRecHits; ++i)
0167 depthValues[i] = rechits[i]->id().depth();
0168 if (nRecHits > 1U)
0169 std::sort(depthValues, depthValues + nRecHits);
0170
0171
0172
0173 uint32_t auxHBHE = 0;
0174 for (unsigned i = 0; i < nRecHits; ++i)
0175 setField(&auxHBHE, 0xf, i * 4, depthValues[i]);
0176 rh->setAuxHBHE(auxHBHE);
0177 }