Initial implementation of Plan 1 rechit combination

RecoLocalCalo/HcalRecAlgos/interface/AbsPlan1RechitCombiner.h

@@ -0,0 +1,42 @@
+#ifndef RecoLocalCalo_HcalRecAlgos_AbsPlan1RechitCombiner_h_
+#define RecoLocalCalo_HcalRecAlgos_AbsPlan1RechitCombiner_h_
+
+#include <utility>
+
+#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
+
+class HcalTopology;
+
+class AbsPlan1RechitCombiner
+{
+public:
+    inline virtual ~AbsPlan1RechitCombiner() {}
+
+    // The topology should be set before the first call
+    // to "add" and whenever it changes
+    virtual void setTopo(const HcalTopology* topo) = 0;
+
+    // The "clear" function is called once per event,
+    // at the beginning of the rechit processing
+    virtual void clear() = 0;
+
+    // This method should be called to add a rechit to process.
+    // It will be assumed that the rechit reference will remain
+    // valid at the time "combine" method is called.
+    virtual void add(const HBHERecHit& rh) = 0;
+
+    // This method should be called once per event,
+    // after all rechits have been added
+    virtual void combine(HBHERecHitCollection* collectionToFill) = 0;
+
+protected:
+    // The first element of the pair is the value to average
+    // and the second is the weight (energy). Non-positive weights
+    // will be ignored.
+    typedef std::pair<float, float> FPair;
+
+    static float energyWeightedAverage(const FPair* data, unsigned len,
+                                       float valueToReturnOnFailure);
+};
+
+#endif // RecoLocalCalo_HcalRecAlgos_AbsPlan1RechitCombiner_h_

RecoLocalCalo/HcalRecAlgos/interface/CaloRecHitAuxSetter.h

@@ -17,6 +17,12 @@ namespace CaloRecHitAuxSetter
     inline void setBit(uint32_t* u, const unsigned bitnum, const bool b)
         {if (b) {*u |= (1U << bitnum);} else {*u &= ~(1U << bitnum);}}
 
+    inline void orBit(uint32_t* u, const unsigned bitnum, const bool b)
+        {if (b) {*u |= (1U << bitnum);}}
+
+    inline void andBit(uint32_t* u, const unsigned bitnum, const bool b)
+        {if (!b) {*u &= ~(1U << bitnum);}}
+
     inline bool getBit(const uint32_t u, const unsigned bitnum)
         {return u & (1U << bitnum);}
 }

RecoLocalCalo/HcalRecAlgos/interface/HBHEIsolatedNoiseAlgos.h

@@ -117,8 +117,9 @@ class ObjectValidator : public ObjectValidatorAbs
 
   uint32_t HcalAcceptSeverityLevel_; // severity level to accept HCAL hits
   uint32_t EcalAcceptSeverityLevel_; // severity level to accept ECAL hits
-  bool UseHcalRecoveredHits_; // whether or not to use recovered HCAL hits
-  bool UseEcalRecoveredHits_; // whether or not to use recovered HCAL hits
+  bool UseHcalRecoveredHits_;  // whether or not to use recovered HCAL hits
+  bool UseEcalRecoveredHits_;  // whether or not to use recovered HCAL hits
+  bool UseAllCombinedRechits_; // whether to use all "Plan 1" combined rechits
 
   double MinValidTrackPt_; // minimum valid track pT
   double MinValidTrackPtBarrel_; // minimum valid track pT in the Barrel

RecoLocalCalo/HcalRecAlgos/interface/HBHERecHitAuxSetter.h

@@ -37,6 +37,9 @@ struct HBHERecHitAuxSetter
     static const unsigned OFF_LINK_ERR = 28;
     static const unsigned OFF_CAPID_ERR = 29;
 
+    // Flag identifying combined "Plan 1" rechits
+    static const unsigned OFF_COMBINED = 30;
+
     // Main function for setting the aux words.
     static void setAux(const HBHEChannelInfo& info,
                        HBHERecHit* rechit);

RecoLocalCalo/HcalRecAlgos/interface/SimplePlan1RechitCombiner.h

@@ -0,0 +1,42 @@
+#ifndef RecoLocalCalo_HcalRecAlgos_SimplePlan1RechitCombiner_h_
+#define RecoLocalCalo_HcalRecAlgos_SimplePlan1RechitCombiner_h_
+
+// Base class header
+#include "RecoLocalCalo/HcalRecAlgos/interface/AbsPlan1RechitCombiner.h"
+
+class SimplePlan1RechitCombiner : public AbsPlan1RechitCombiner
+{
+public:
+    SimplePlan1RechitCombiner();
+
+    inline virtual ~SimplePlan1RechitCombiner() {}
+
+    virtual void setTopo(const HcalTopology* topo) override;
+    virtual void clear() override;
+    virtual void add(const HBHERecHit& rh) override;
+    virtual void combine(HBHERecHitCollection* toFill) override;
+
+protected:
+    // Map the original detector id into the id of the composite rechit
+    virtual HcalDetId mapRechit(HcalDetId from) const;
+
+    // Return rechit with id of 0 if it is to be discarded
+    virtual HBHERecHit makeRechit(HcalDetId idToMake,
+                                  const std::vector<const HBHERecHit*>& rechits) const;
+
+    // Combine rechit auxiliary information
+    virtual void combineAuxInfo(const std::vector<const HBHERecHit*>& rechits,
+                                HBHERecHit* rh) const;
+
+    const HcalTopology* topo_;
+
+private:
+    // The first element of the pair is the id of the combined rechit.
+    // The second element is the pointer to the original rechit.
+    typedef std::pair<HcalDetId, const HBHERecHit*> MapItem;
+
+    std::vector<MapItem> rechitMap_;
+    std::vector<const HBHERecHit*> ptrbuf_;
+};
+
+#endif // RecoLocalCalo_HcalRecAlgos_SimplePlan1RechitCombiner_h_

RecoLocalCalo/HcalRecAlgos/interface/parsePlan1RechitCombiner.h

@@ -0,0 +1,18 @@
+#ifndef RecoLocalCalo_HcalRecAlgos_parsePlan1RechitCombiner_h_
+#define RecoLocalCalo_HcalRecAlgos_parsePlan1RechitCombiner_h_
+
+#include <memory>
+#include "RecoLocalCalo/HcalRecAlgos/interface/AbsPlan1RechitCombiner.h"
+
+namespace edm {
+    class ParameterSet;
+}
+
+//
+// Factory function for creating objects of types inheriting from
+// AbsPlan1RechitCombiner out of parameter sets
+//
+std::unique_ptr<AbsPlan1RechitCombiner>
+parsePlan1RechitCombiner(const edm::ParameterSet& ps);
+
+#endif // RecoLocalCalo_HcalRecAlgos_parsePlan1RechitCombiner_h_

RecoLocalCalo/HcalRecAlgos/src/AbsPlan1RechitCombiner.cc

@@ -0,0 +1,21 @@
+#include "RecoLocalCalo/HcalRecAlgos/interface/AbsPlan1RechitCombiner.h"
+
+float AbsPlan1RechitCombiner::energyWeightedAverage(
+    const FPair* data, const unsigned len,
+    const float valueToReturnOnFailure)
+{
+    double sum = 0.0, wsum = 0.0;
+    for (unsigned i=0; i<len; ++i)
+    {
+        const float w = data[i].second;
+        if (w > 0.f)
+        {
+            sum += w*data[i].first;
+            wsum += w;
+        }
+    }
+    if (wsum > 0.0)
+        return sum/wsum;
+    else
+        return valueToReturnOnFailure;
+}

RecoLocalCalo/HcalRecAlgos/src/HBHEIsolatedNoiseAlgos.cc

@@ -31,6 +31,8 @@ Original Author: John Paul Chou (Brown University)
 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
+#include "RecoLocalCalo/HcalRecAlgos/interface/HBHERecHitAuxSetter.h"
+#include "RecoLocalCalo/HcalRecAlgos/interface/CaloRecHitAuxSetter.h"
 
 ////////////////////////////////////////////////////////////
 //
@@ -50,6 +52,7 @@ ObjectValidator::ObjectValidator(const edm::ParameterSet& iConfig)
   EcalAcceptSeverityLevel_ = iConfig.getParameter<uint32_t>("EcalAcceptSeverityLevel");
   UseHcalRecoveredHits_ = iConfig.getParameter<bool>("UseHcalRecoveredHits");
   UseEcalRecoveredHits_ = iConfig.getParameter<bool>("UseEcalRecoveredHits");
+  UseAllCombinedRechits_ = iConfig.getParameter<bool>("UseAllCombinedRechits");
 
   MinValidTrackPt_ = iConfig.getParameter<double>("MinValidTrackPt");
   MinValidTrackPtBarrel_ = iConfig.getParameter<double>("MinValidTrackPtBarrel");
@@ -73,6 +76,10 @@ bool ObjectValidator::validHit(const HBHERecHit& hit) const
 {
   assert(theHcalSevLvlComputer_!=0 && theHcalChStatus_!=0);
 
+  if (UseAllCombinedRechits_)
+      if (CaloRecHitAuxSetter::getBit(hit.auxPhase1(), HBHERecHitAuxSetter::OFF_COMBINED))
+          return true;
+
   // require the hit to pass a certain energy threshold
   if(hit.id().subdet()==HcalBarrel && hit.energy()<HBThreshold_) return false;
   else if(hit.id().subdet()==HcalEndcap && hit.id().ietaAbs()<=20 && hit.energy()<HESThreshold_) return false;

RecoLocalCalo/HcalRecAlgos/src/SimplePlan1RechitCombiner.cc

@@ -0,0 +1,188 @@
+#include <cassert>
+#include <algorithm>
+
+#include "DataFormats/HcalRecHit/interface/HcalSpecialTimes.h"
+
+#include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"
+
+#include "Geometry/CaloTopology/interface/HcalTopology.h"
+
+#include "RecoLocalCalo/HcalRecAlgos/interface/SimplePlan1RechitCombiner.h"
+#include "RecoLocalCalo/HcalRecAlgos/interface/HBHERecHitAuxSetter.h"
+#include "RecoLocalCalo/HcalRecAlgos/interface/CaloRecHitAuxSetter.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;
+
+    const unsigned nRecHits = rechits.size();
+    assert(nRecHits);
+    assert(nRecHits <= HBHERecHitAuxSetter::MASK_NSAMPLES);
+
+    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);
+}

RecoLocalCalo/HcalRecAlgos/src/parsePlan1RechitCombiner.cc

@@ -0,0 +1,18 @@
+#include "RecoLocalCalo/HcalRecAlgos/interface/parsePlan1RechitCombiner.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+// Concrete "Plan 1" rechit combination algorithm headers
+#include "RecoLocalCalo/HcalRecAlgos/interface/SimplePlan1RechitCombiner.h"
+
+std::unique_ptr<AbsPlan1RechitCombiner>
+parsePlan1RechitCombiner(const edm::ParameterSet& ps)
+{
+    std::unique_ptr<AbsPlan1RechitCombiner> algo;
+
+    const std::string& className = ps.getParameter<std::string>("Class");
+
+    if (className == "SimplePlan1RechitCombiner")
+        algo = std::unique_ptr<AbsPlan1RechitCombiner>(new SimplePlan1RechitCombiner());
+
+    return algo;
+}