HcalRecAlgos/src/SimplePlan1RechitCombiner.cc

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;  // Should be >= max # of Phase 1 HCAL depths
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   // Combine energies, times, and fit chi-square
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   // Combine the auxiliary information
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   // The number of rechits should be not larger than the
0104   // number of half-bytes in a 32-bit word
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   // Combine various status bits
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     // Status flags are simply ORed for now. Might want
0135     // to rethink this in the future.
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   // A number that can be later used to calculate chi-square NDoF
0148   setField(&auxPhase1, HBHERecHitAuxSetter::MASK_NSAMPLES, HBHERecHitAuxSetter::OFF_ADC, tripleFitCount);
0149
0150   // How many rechits were combined?
0151   setField(&auxPhase1, HBHERecHitAuxSetter::MASK_NSAMPLES, HBHERecHitAuxSetter::OFF_NSAMPLES, nRecHits);
0152
0153   // Should combine QIE11 data only
0154   setBit(&auxPhase1, HBHERecHitAuxSetter::OFF_TDC_TIME, true);
0155
0156   // Indicate that this rechit is combined
0157   setBit(&auxPhase1, HBHERecHitAuxSetter::OFF_COMBINED, true);
0158
0159   // Copy the aux words into the rechit
0160   rh->setFlags(flags);
0161   rh->setAuxPhase1(auxPhase1);
0162
0163   // Sort the depth values of the combined rechits
0164   // in the increasing order
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   // Pack the information about the depth of the rechits
0172   // that we are combining into the "auxHBHE" word
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 }