SimplePlan1RechitCombiner.cc

Go to the documentation of this file.

#include <cassert>
#include <algorithm>

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

#include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"

#include "Geometry/CaloTopology/interface/HcalTopology.h"

#include "RecoLocalCalo/HcalRecAlgos/interface/SimplePlan1RechitCombiner.h"

SimplePlan1RechitCombiner::SimplePlan1RechitCombiner() : topo_(nullptr) {}

void SimplePlan1RechitCombiner::setTopo(const HcalTopology* topo) { topo_ = topo; }

void SimplePlan1RechitCombiner::clear() { rechitMap_.clear(); }

void SimplePlan1RechitCombiner::add(const HBHERecHit& rh) {
  if (!CaloRecHitAuxSetter::getBit(rh.auxPhase1(), HBHERecHitAuxSetter::OFF_DROPPED))
    rechitMap_.push_back(MapItem(mapRechit(rh.id()), &rh));
}

void SimplePlan1RechitCombiner::combine(HBHERecHitCollection* toFill) {
  if (!rechitMap_.empty()) {
    std::sort(rechitMap_.begin(), rechitMap_.end());

    HcalDetId oldId(rechitMap_[0].first);
    ptrbuf_.clear();
    ptrbuf_.push_back(rechitMap_[0].second);

    const std::size_t nInput = rechitMap_.size();
    for (std::size_t i = 1; i < nInput; ++i) {
      if (rechitMap_[i].first != oldId) {
        const HBHERecHit& rh = makeRechit(oldId, ptrbuf_);
        if (rh.id().rawId())
          toFill->push_back(rh);
        oldId = rechitMap_[i].first;
        ptrbuf_.clear();
      }
      ptrbuf_.push_back(rechitMap_[i].second);
    }

    const HBHERecHit& rh = makeRechit(oldId, ptrbuf_);
    if (rh.id().rawId())
      toFill->push_back(rh);
  }
}

HcalDetId SimplePlan1RechitCombiner::mapRechit(HcalDetId from) const { return topo_->mergedDepthDetId(from); }

HBHERecHit SimplePlan1RechitCombiner::makeRechit(const HcalDetId idToMake,
                                                 const std::vector<const HBHERecHit*>& rechits) const {
  constexpr unsigned MAXLEN = 8U;  // Should be >= max # of Phase 1 HCAL depths
  constexpr float TIME_IF_NO_ENERGY = -999.f;

  const unsigned nRecHits = rechits.size();
  assert(nRecHits);
  assert(nRecHits <= MAXLEN);

  // Combine energies, times, and fit chi-square
  double energy = 0.0, eraw = 0.0, eaux = 0.0, chisq = 0.0;
  FPair times[MAXLEN], adctimes[MAXLEN];
  unsigned nADCTimes = 0;

  for (unsigned i = 0; i < nRecHits; ++i) {
    const HBHERecHit& rh(*rechits[i]);
    const float e = rh.energy();
    energy += e;
    eraw += rh.eraw();
    eaux += rh.eaux();
    chisq += rh.chi2();
    times[i].first = rh.time();
    times[i].second = e;

    const float tADC = rh.timeFalling();
    if (!HcalSpecialTimes::isSpecial(tADC)) {
      adctimes[nADCTimes].first = tADC;
      adctimes[nADCTimes].second = e;
      ++nADCTimes;
    }
  }

  HBHERecHit rh(idToMake,
                energy,
                energyWeightedAverage(times, nRecHits, TIME_IF_NO_ENERGY),
                energyWeightedAverage(adctimes, nADCTimes, HcalSpecialTimes::UNKNOWN_T_NOTDC));
  rh.setRawEnergy(eraw);
  rh.setAuxEnergy(eaux);
  rh.setChiSquared(chisq);

  // Combine the auxiliary information
  combineAuxInfo(rechits, &rh);

  return rh;
}

void SimplePlan1RechitCombiner::combineAuxInfo(const std::vector<const HBHERecHit*>& rechits, HBHERecHit* rh) const {
  using namespace CaloRecHitAuxSetter;
  using namespace HcalPhase1FlagLabels;

  // The number of rechits should be not larger than the
  // number of half-bytes in a 32-bit word
  constexpr unsigned MAXLEN = 8U;

  const unsigned nRecHits = rechits.size();
  assert(nRecHits);
  assert(nRecHits <= MAXLEN);

  uint32_t flags = 0, auxPhase1 = 0;
  unsigned tripleFitCount = 0;
  unsigned soiVote[HBHERecHitAuxSetter::MASK_SOI + 1U] = {0};
  unsigned capidVote[HBHERecHitAuxSetter::MASK_CAPID + 1U] = {0};

  // Combine various status bits
  for (unsigned i = 0; i < nRecHits; ++i) {
    const HBHERecHit& rh(*rechits[i]);
    const uint32_t rhflags = rh.flags();
    const uint32_t rhAuxPhase1 = rh.auxPhase1();

    orBit(&auxPhase1, HBHERecHitAuxSetter::OFF_LINK_ERR, getBit(rhAuxPhase1, HBHERecHitAuxSetter::OFF_LINK_ERR));
    orBit(&auxPhase1, HBHERecHitAuxSetter::OFF_CAPID_ERR, getBit(rhAuxPhase1, HBHERecHitAuxSetter::OFF_CAPID_ERR));

    const unsigned soi = getField(rhAuxPhase1, HBHERecHitAuxSetter::MASK_SOI, HBHERecHitAuxSetter::OFF_SOI);
    soiVote[soi]++;

    const unsigned capid = getField(rhAuxPhase1, HBHERecHitAuxSetter::MASK_CAPID, HBHERecHitAuxSetter::OFF_CAPID);
    capidVote[capid]++;

    if (getBit(rhflags, HBHEStatusFlag::HBHEPulseFitBit))
      ++tripleFitCount;

    // Status flags are simply ORed for now. Might want
    // to rethink this in the future.
    flags |= rhflags;
  }

  unsigned* pmaxsoi = std::max_element(soiVote, soiVote + sizeof(soiVote) / sizeof(soiVote[0]));
  const unsigned soi = std::distance(&soiVote[0], pmaxsoi);
  setField(&auxPhase1, HBHERecHitAuxSetter::MASK_SOI, HBHERecHitAuxSetter::OFF_SOI, soi);

  unsigned* pmaxcapid = std::max_element(capidVote, capidVote + sizeof(capidVote) / sizeof(capidVote[0]));
  const unsigned capid = std::distance(&capidVote[0], pmaxcapid);
  setField(&auxPhase1, HBHERecHitAuxSetter::MASK_CAPID, HBHERecHitAuxSetter::OFF_CAPID, capid);

  // A number that can be later used to calculate chi-square NDoF
  setField(&auxPhase1, HBHERecHitAuxSetter::MASK_NSAMPLES, HBHERecHitAuxSetter::OFF_ADC, tripleFitCount);

  // How many rechits were combined?
  setField(&auxPhase1, HBHERecHitAuxSetter::MASK_NSAMPLES, HBHERecHitAuxSetter::OFF_NSAMPLES, nRecHits);

  // Should combine QIE11 data only
  setBit(&auxPhase1, HBHERecHitAuxSetter::OFF_TDC_TIME, true);

  // Indicate that this rechit is combined
  setBit(&auxPhase1, HBHERecHitAuxSetter::OFF_COMBINED, true);

  // Copy the aux words into the rechit
  rh->setFlags(flags);
  rh->setAuxPhase1(auxPhase1);

  // Sort the depth values of the combined rechits
  // in the increasing order
  unsigned depthValues[MAXLEN];
  for (unsigned i = 0; i < nRecHits; ++i)
    depthValues[i] = rechits[i]->id().depth();
  if (nRecHits > 1U)
    std::sort(depthValues, depthValues + nRecHits);

  // Pack the information about the depth of the rechits
  // that we are combining into the "auxHBHE" word
  uint32_t auxHBHE = 0;
  for (unsigned i = 0; i < nRecHits; ++i)
    setField(&auxHBHE, 0xf, i * 4, depthValues[i]);
  rh->setAuxHBHE(auxHBHE);
}