Infrastructure for using channel-dependent pulse shapes for HB/HE

CalibCalorimetry/HcalAlgos/interface/HcalPulseShapeLookup.h

@@ -0,0 +1,43 @@
+#ifndef CalibCalorimetry_HcalAlgos_HcalPulseShapeLookup_h
+#define CalibCalorimetry_HcalAlgos_HcalPulseShapeLookup_h
+
+#include <vector>
+#include <string>
+#include <utility>
+
+#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShape.h"
+#include "Geometry/CaloTopology/interface/HcalTopology.h"
+
+// This class must be sufficiently similar in its interface to HcalPulseShapes
+// so that both of them can be used as parameters of the same templated code.
+// However, this one is designed for use with a more efficient pulse shape
+// lookup scheme and can accommodate more pulse shapes.
+class HcalPulseShapeLookup {
+public:
+  typedef HcalPulseShape Shape;
+  typedef std::pair<std::string, Shape> LabeledShape;
+
+  HcalPulseShapeLookup(const std::vector<LabeledShape>& shapes,
+                       const std::vector<int>& channelToTypeLookup,
+                       const HcalTopology* htopo);
+
+  inline unsigned nShapeTypes() const { return theShapes_.size(); }
+  const Shape& getShape(int shapeType) const;
+  const std::string& getLabel(int shapeType) const;
+
+  int getShapeType(unsigned linearizedChannelNumber) const;
+  const Shape& getChannelShape(unsigned linearizedChannelNumber) const;
+  const std::string& getChannelLabel(unsigned linearizedChannelNumber) const;
+
+  int getShapeType(const DetId& id) const;
+  const Shape& getChannelShape(const DetId& id) const;
+  const std::string& getChannelLabel(const DetId& id) const;
+
+private:
+  std::vector<LabeledShape> theShapes_;
+  std::vector<int> shapeTypes_;
+  // We do not own the pointer
+  const HcalTopology* htopo_;
+};
+
+#endif  // CalibCalorimetry_HcalAlgos_HcalPulseShapeLookup_h

CalibCalorimetry/HcalAlgos/src/ES_HcalPulseShapeLookup.cc

@@ -0,0 +1,4 @@
+#include "FWCore/Utilities/interface/typelookup.h"
+#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapeLookup.h"
+
+TYPELOOKUP_DATA_REG(HcalPulseShapeLookup);

CalibCalorimetry/HcalAlgos/src/HcalPulseShapeLookup.cc

@@ -0,0 +1,38 @@
+#include <cassert>
+
+#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapeLookup.h"
+
+HcalPulseShapeLookup::HcalPulseShapeLookup(const std::vector<LabeledShape>& shapes,
+                                           const std::vector<int>& channelToTypeLookup,
+                                           const HcalTopology* htopo)
+    : theShapes_(shapes), shapeTypes_(channelToTypeLookup), htopo_(htopo) {
+  assert(htopo_);
+}
+
+const HcalPulseShapeLookup::Shape& HcalPulseShapeLookup::getShape(const int shapeType) const {
+  return theShapes_.at(shapeType).second;
+}
+
+const std::string& HcalPulseShapeLookup::getLabel(const int shapeType) const { return theShapes_.at(shapeType).first; }
+
+int HcalPulseShapeLookup::getShapeType(const unsigned linearizedChannelNumber) const {
+  return shapeTypes_.at(linearizedChannelNumber);
+}
+
+int HcalPulseShapeLookup::getShapeType(const DetId& id) const { return getShapeType(htopo_->detId2denseId(id)); }
+
+const HcalPulseShapeLookup::Shape& HcalPulseShapeLookup::getChannelShape(const unsigned linearizedChannelNumber) const {
+  return theShapes_.at(shapeTypes_.at(linearizedChannelNumber)).second;
+}
+
+const std::string& HcalPulseShapeLookup::getChannelLabel(const unsigned linearizedChannelNumber) const {
+  return theShapes_.at(shapeTypes_.at(linearizedChannelNumber)).first;
+}
+
+const HcalPulseShapeLookup::Shape& HcalPulseShapeLookup::getChannelShape(const DetId& id) const {
+  return getChannelShape(htopo_->detId2denseId(id));
+}
+
+const std::string& HcalPulseShapeLookup::getChannelLabel(const DetId& id) const {
+  return getChannelLabel(htopo_->detId2denseId(id));
+}

DataFormats/HcalRecHit/test/HcalRecHitDump.cc

@@ -87,6 +87,9 @@ namespace {
     } else {
       s << " none";
     }
+
+    // Dump fit chi^2. It will be set to -1.0f if it is not available.
+    s << "; chi2: " << i.chi2();
   }
 
   void printRecHitAuxInfo(std::ostream& s, const HFRecHit& i, const std::vector<int>& bits, bool) {

RecoLocalCalo/Configuration/python/hcalLocalReco_cff.py

@@ -3,6 +3,7 @@
 
 from RecoLocalCalo.HcalRecAlgos.hcalRecAlgoESProd_cfi import *
 from RecoLocalCalo.HcalRecAlgos.hcalChannelPropertiesESProd_cfi import *
+from RecoLocalCalo.HcalRecAlgos.hcalPulseShapesESProd_cfi import *
 hcalOOTPileupESProducer = cms.ESProducer('OOTPileupDBCompatibilityESProducer')
 
 from RecoLocalCalo.HcalRecProducers.HBHEPhase1Reconstructor_cfi import hbheprereco as _phase1_hbheprereco

RecoLocalCalo/HcalRecAlgos/interface/AbsHBHEPhase1Algo.h

@@ -5,7 +5,9 @@
 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
 #include "DataFormats/HcalRecHit/interface/HBHEChannelInfo.h"
 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
+#include "Geometry/CaloTopology/interface/HcalTopology.h"
 #include "CondFormats/HcalObjects/interface/HcalRecoParam.h"
+#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapeLookup.h"
 
 class AbsHcalAlgoData;
 
@@ -26,6 +28,7 @@ class AbsHcalAlgoData;
 //
 class AbsHBHEPhase1Algo {
 public:
+  inline AbsHBHEPhase1Algo() : channelPulseShapes_(nullptr) {}
   inline virtual ~AbsHBHEPhase1Algo() {}
 
   inline virtual void beginRun(const edm::Run&, const edm::EventSetup&) {}
@@ -50,6 +53,11 @@ class AbsHBHEPhase1Algo {
                                  const HcalRecoParam* params,
                                  const HcalCalibrations& calibs,
                                  bool isRealData) = 0;
+
+  inline virtual void setPulseShapes(const HcalPulseShapeLookup* ptr) { channelPulseShapes_ = ptr; }
+
+protected:
+  const HcalPulseShapeLookup* channelPulseShapes_;
 };
 
 #endif  // RecoLocalCalo_HcalRecAlgos_AbsHBHEPhase1Algo_h_

RecoLocalCalo/HcalRecAlgos/interface/HcalPulseShapeLookupRcd.h

@@ -0,0 +1,12 @@
+#ifndef RecoLocalCalo_HcalRecAlgos_HcalPulseShapeLookupRcd_h
+#define RecoLocalCalo_HcalRecAlgos_HcalPulseShapeLookupRcd_h
+
+#include "FWCore/Framework/interface/DependentRecordImplementation.h"
+#include "Geometry/Records/interface/HcalRecNumberingRecord.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+
+class HcalPulseShapeLookupRcd : public edm::eventsetup::DependentRecordImplementation<
+                                    HcalPulseShapeLookupRcd,
+                                    edm::mpl::Vector<HcalRecNumberingRecord, CaloGeometryRecord> > {};
+
+#endif  // RecoLocalCalo_HcalRecAlgos_HcalPulseShapeLookupRcd_h

RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h

@@ -96,7 +96,7 @@ struct MahiDebugInfo {
 class MahiFit {
 public:
   MahiFit();
-  ~MahiFit() {}
+  ~MahiFit() {};
 
   void setParameters(bool iDynamicPed,
                      double iTS4Thresh,
@@ -126,12 +126,38 @@ class MahiFit {
 
   void doFit(std::array<float, 4>& correctedOutput, const int nbx) const;
 
-  void setPulseShapeTemplate(int pulseShapeId,
-                             const HcalPulseShapes& ps,
-                             bool hasTimeInfo,
+  // PulseShapeCollection is expected to be either HcalPulseShapes
+  // or HcalPulseShapeLookup
+  template <class PulseShapeCollection>
+  void setPulseShapeTemplate(const int pulseShapeId,
+                             const PulseShapeCollection& ps,
+                             const bool hasTimeInfo,
                              const HcalTimeSlew* hcalTimeSlewDelay,
-                             unsigned int nSamples,
-                             const float gain);
+                             const unsigned int nSamples,
+                             const float gain0) {
+    if (hcalTimeSlewDelay != hcalTimeSlewDelay_) {
+      assert(hcalTimeSlewDelay);
+      hcalTimeSlewDelay_ = hcalTimeSlewDelay;
+      tsDelay1GeV_ = hcalTimeSlewDelay->delay(1.0, slewFlavor_);
+    }
+
+    if (pulseShapeId != currentPulseShapeId_) {
+      resetPulseShapeTemplate(pulseShapeId, ps, nSamples);
+    }
+
+    // 1 sigma time constraint
+    nnlsWork_.dt = hasTimeInfo ? timeSigmaSiPM_ : timeSigmaHPD_;
+
+    if (nnlsWork_.tsSize != nSamples) {
+      nnlsWork_.tsSize = nSamples;
+      nnlsWork_.amplitudes.resize(nSamples);
+      nnlsWork_.noiseTerms.resize(nSamples);
+      nnlsWork_.pedVals.resize(nSamples);
+    }
+
+    // threshold in GeV for ccTime
+    thEnergeticPulsesFC_ = thEnergeticPulses_ / gain0;
+  }
 
   typedef BXVector::Index Index;
   const HcalTimeSlew* hcalTimeSlewDelay_ = nullptr;
@@ -145,7 +171,30 @@ class MahiFit {
   const float minimize() const;
   void onePulseMinimize() const;
   void updateCov(const SampleMatrix& invCovMat) const;
-  void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples);
+
+  template <class PulseShapeCollection>
+  void resetPulseShapeTemplate(const int pulseShapeId, const PulseShapeCollection& ps, unsigned int /* nSamples */) {
+    ++cntsetPulseShape_;
+
+    psfPtr_ = nullptr;
+    for (auto& elem : knownPulseShapes_) {
+      if (elem.first == pulseShapeId) {
+        psfPtr_ = &*elem.second;
+        break;
+      }
+    }
+
+    if (!psfPtr_) {
+      // only the pulse shape itself from PulseShapeFunctor is used for Mahi
+      // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
+      auto p = std::make_shared<FitterFuncs::PulseShapeFunctor>(
+          ps.getShape(pulseShapeId), false, false, false, 1, 0, 0, hcal::constants::maxSamples);
+      knownPulseShapes_.emplace_back(pulseShapeId, p);
+      psfPtr_ = &*p;
+    }
+
+    currentPulseShapeId_ = pulseShapeId;
+  }
 
   float ccTime(const float itQ) const;
   void updatePulseShape(const float itQ,