Infrastructure for HCAL channel-dependent pulse shapes using database configurations

CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h

@@ -22,6 +22,8 @@ Text file formats for different data types is as following:
   eta(int)  phi(int) depth(int) det(HB,HE,HF) cap1_value(float) cap2_value(float) cap3_value(float) cap4_value(float)  HcalDetId(int,optional)
 - HcalPFCuts:
   eta(int)  phi(int) depth(int) det(HB,HE,HF) noiseThreshold(float) seedThreshold(float)
+- HcalPulseDelays:
+  eta(int)  phi(int) depth(int) det(HB,HE,HF) delay(float)
 - HcalPedestalWidths:
   eta(int)  phi(int) depth(int) det(HB,HE,HF) sigma_1_1(float) sigma_2_1 sigma_2_2 sigma_3_1 sigma_3_2 sigma_3_3 sigma_4_1 sigma_4_2 sigma_4_3 sigma_4_4
 - HcalQIEShape:
@@ -71,6 +73,8 @@ namespace HcalDbASCIIIO {
   bool dumpObject(std::ostream& fOutput, const HcalGainWidths& fObject);
   bool getObject(std::istream& fInput, HcalPFCuts* fObject);
   bool dumpObject(std::ostream& fOutput, const HcalPFCuts& fObject);
+  bool getObject(std::istream& fInput, HcalPulseDelays* fObject);
+  bool dumpObject(std::ostream& fOutput, const HcalPulseDelays& fObject);
   bool getObject(std::istream& fInput, HcalQIEData* fObject);
   bool dumpObject(std::ostream& fOutput, const HcalQIEData& fObject);
   bool getObject(std::istream& fInput, HcalCalibrationQIEData* fObject);

CalibCalorimetry/HcalAlgos/interface/HcalDbHardcode.h

@@ -12,6 +12,7 @@
 #include "CondFormats/HcalObjects/interface/HcalGain.h"
 #include "CondFormats/HcalObjects/interface/HcalGainWidth.h"
 #include "CondFormats/HcalObjects/interface/HcalPFCut.h"
+#include "CondFormats/HcalObjects/interface/HcalPulseDelay.h"
 #include "CondFormats/HcalObjects/interface/HcalZSThreshold.h"
 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
 #include "CondFormats/HcalObjects/interface/HcalQIEShape.h"
@@ -102,6 +103,7 @@ class HcalDbHardcode {
   HcalGain makeGain(HcalGenericDetId fId, bool fSmear = false) const;
   HcalGainWidth makeGainWidth(HcalGenericDetId fId) const;
   HcalPFCut makePFCut(HcalGenericDetId fId, double intlumi, bool noHE) const;
+  HcalPulseDelay makePulseDelay(HcalGenericDetId fId) const;
   HcalZSThreshold makeZSThreshold(HcalGenericDetId fId) const;
   HcalQIECoder makeQIECoder(HcalGenericDetId fId) const;
   HcalCalibrationQIECoder makeCalibrationQIECoder(HcalGenericDetId fId) const;

CalibCalorimetry/HcalAlgos/interface/HcalHardcodeParameters.h

@@ -26,7 +26,8 @@ class HcalHardcodeParameters {
                          const std::vector<double>& darkCurrent,
                          const std::vector<double>& noiseCorrelation,
                          double noiseThreshold,
-                         double seedThreshold);
+                         double seedThreshold,
+                         double pulseDelay);
 
   //construct from pset
   HcalHardcodeParameters(const edm::ParameterSet& p);
@@ -51,6 +52,7 @@ class HcalHardcodeParameters {
   double noiseCorrelation(unsigned index) const;
   inline double noiseThreshold() const { return noiseThreshold_; }
   inline double seedThreshold() const { return seedThreshold_; }
+  inline double pulseDelay() const { return pulseDelay_; }
 
 private:
   //member variables
@@ -66,6 +68,7 @@ class HcalHardcodeParameters {
   bool doSipmRadiationDamage_;
   HcalSiPMRadiationDamage sipmRadiationDamage_;
   double noiseThreshold_, seedThreshold_;
+  double pulseDelay_;
 };
 
 #endif

CalibCalorimetry/HcalAlgos/interface/HcalPulseShapeLookup.h

@@ -0,0 +1,48 @@
+#ifndef CalibCalorimetry_HcalAlgos_HcalPulseShapeLookup_h
+#define CalibCalorimetry_HcalAlgos_HcalPulseShapeLookup_h
+
+#include <vector>
+#include <string>
+#include <tuple>
+
+#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::tuple<std::string, float, 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;
+  float getTimeShift(int shapeType) const;
+
+  int getShapeType(unsigned linearizedChannelNumber) const;
+  const Shape& getChannelShape(unsigned linearizedChannelNumber) const;
+  const std::string& getChannelLabel(unsigned linearizedChannelNumber) const;
+  float getChannelTimeShift(unsigned linearizedChannelNumber) const;
+
+  int getShapeType(const DetId& id) const;
+  const Shape& getChannelShape(const DetId& id) const;
+  const std::string& getChannelLabel(const DetId& id) const;
+  float getChannelTimeShift(const DetId& id) const;
+
+  void dumpToTxt(const std::string& filename, unsigned precision = 0U) 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/HcalDbASCIIIO.cc

@@ -56,8 +56,9 @@ std::vector<std::string> splitString(const std::string& fLine) {
   std::vector<std::string> result;
   int start = 0;
   bool empty = true;
-  for (unsigned i = 0; i <= fLine.size(); i++) {
-    if (fLine[i] == ' ' || i == fLine.size()) {
+  const unsigned sz = fLine.size();
+  for (unsigned i = 0; i <= sz; i++) {
+    if (i == sz || fLine[i] == ' ' || fLine[i] == '\t') {
       if (!empty) {
         std::string item(fLine, start, i - start);
         result.push_back(item);
@@ -303,6 +304,49 @@ bool getHcalDoubleFloatObject(std::istream& fInput, T* fObject) {
   return true;
 }
 
+template <class T>
+bool dumpHcalStringFloatObject(std::ostream& fOutput, const T& fObject) {
+  char buffer[1024];
+  sprintf(buffer, "# %15s %15s %15s %15s %10s %10s %10s\n", "eta", "phi", "dep", "det", "value0", "value1", "DetId");
+  fOutput << buffer;
+  std::vector<DetId> channels = fObject.getAllChannels();
+  std::sort(channels.begin(), channels.end(), DetIdLess());
+  for (std::vector<DetId>::iterator channel = channels.begin(); channel != channels.end(); ++channel) {
+    const std::string& value0 = fObject.getValues(*channel)->getValue0();
+    const float value1 = fObject.getValues(*channel)->getValue1();
+    HcalDbASCIIIO::dumpId(fOutput, *channel);
+    fOutput << ' ' << value0;
+    sprintf(buffer, " %10.7f %10X\n", value1, channel->rawId());
+    fOutput << buffer;
+  }
+  return true;
+}
+
+template <class S, class T>
+bool getHcalStringFloatObject(std::istream& fInput, T* fObject) {
+  if (!fObject)
+    return false;  //fObject = new T;
+  char buffer[1024];
+  while (fInput.getline(buffer, 1024)) {
+    if (buffer[0] == '#')
+      continue;  //ignore comment
+    std::vector<std::string> items = splitString(std::string(buffer));
+    if (items.empty())
+      continue;  // blank line
+    if (items.size() < 6) {
+      edm::LogWarning("Format Error") << "Bad line: " << buffer
+                                      << "\n line must contain 6 items: eta, phi, depth, subdet, value0, value1"
+                                      << std::endl;
+      continue;
+    }
+    DetId id = HcalDbASCIIIO::getId(items);
+    S fCondObject(id, items[4], atof(items[5].c_str()));
+    fObject->addValues(fCondObject);
+  }
+
+  return true;
+}
+
 template <class T>
 bool dumpHcalSingleFloatObject(std::ostream& fOutput, const T& fObject) {
   char buffer[1024];
@@ -438,6 +482,13 @@ namespace HcalDbASCIIIO {
     return dumpHcalDoubleFloatObject(fOutput, fObject);
   }
 
+  bool getObject(std::istream& fInput, HcalPulseDelays* fObject) {
+    return getHcalStringFloatObject<HcalPulseDelay>(fInput, fObject);
+  }
+  bool dumpObject(std::ostream& fOutput, const HcalPulseDelays& fObject) {
+    return dumpHcalStringFloatObject(fOutput, fObject);
+  }
+
   bool getObject(std::istream& fInput, HcalRespCorrs* fObject) {
     return getHcalSingleObject<float, HcalRespCorr>(fInput, fObject);
   }

CalibCalorimetry/HcalAlgos/src/HcalDbHardcode.cc

@@ -31,7 +31,8 @@ HcalDbHardcode::HcalDbHardcode()
                             {0.0},                 //dark current
                             {0.0},                 //noise correlation
                             0.0,                   //PF noise threshold
-                            0.1                    //PF seed threshold
+                            0.1,                   //PF seed threshold
+                            0.0                    //Extra pulse delay
                             ),
       setHB_(false),
       setHE_(false),
@@ -178,6 +179,11 @@ HcalGainWidth HcalDbHardcode::makeGainWidth(HcalGenericDetId fId) const {  // Ge
   return result;
 }
 
+HcalPulseDelay HcalDbHardcode::makePulseDelay(HcalGenericDetId fId) const {  // ns
+  const float value0 = getParameters(fId).pulseDelay();
+  return HcalPulseDelay(fId, "default", value0);
+}
+
 HcalPFCut HcalDbHardcode::makePFCut(HcalGenericDetId fId, double intLumi, bool noHE) const {  // GeV
 
   // assign default dummy parameters

CalibCalorimetry/HcalAlgos/src/HcalHardcodeParameters.cc

@@ -14,7 +14,8 @@ HcalHardcodeParameters::HcalHardcodeParameters(const double pedestal,
                                                const std::vector<double>& darkCurrent,
                                                const std::vector<double>& noiseCorrelation,
                                                const double noiseTh,
-                                               const double seedTh)
+                                               const double seedTh,
+                                               const double delay)
     : pedestal_(pedestal),
       pedestalWidth_(pedestalWidth),
       gain_(gain),
@@ -30,7 +31,8 @@ HcalHardcodeParameters::HcalHardcodeParameters(const double pedestal,
       noiseCorrelation_(noiseCorrelation),
       doSipmRadiationDamage_(false),
       noiseThreshold_(noiseTh),
-      seedThreshold_(seedTh) {}
+      seedThreshold_(seedTh),
+      pulseDelay_(delay) {}
 
 HcalHardcodeParameters::HcalHardcodeParameters(const edm::ParameterSet& p)
     : pedestal_(p.getParameter<double>("pedestal")),
@@ -48,7 +50,8 @@ HcalHardcodeParameters::HcalHardcodeParameters(const edm::ParameterSet& p)
       noiseCorrelation_(p.getParameter<std::vector<double>>("noiseCorrelation")),
       doSipmRadiationDamage_(p.getParameter<bool>("doRadiationDamage")),
       noiseThreshold_(p.getParameter<double>("noiseThreshold")),
-      seedThreshold_(p.getParameter<double>("seedThreshold")) {
+      seedThreshold_(p.getParameter<double>("seedThreshold")),
+      pulseDelay_(p.getParameter<double>("pulseDelay")) {
   if (doSipmRadiationDamage_)
     sipmRadiationDamage_ = HcalSiPMRadiationDamage(darkCurrent_, p.getParameter<edm::ParameterSet>("radiationDamage"));
 }

CalibCalorimetry/HcalAlgos/src/HcalPulseShapeLookup.cc

@@ -0,0 +1,97 @@
+#include <cassert>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <algorithm>
+
+#include "FWCore/Utilities/interface/Exception.h"
+
+#include "DataFormats/HcalDetId/interface/HcalDetId.h"
+#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 std::get<2>(theShapes_.at(shapeType));
+}
+
+const std::string& HcalPulseShapeLookup::getLabel(const int shapeType) const {
+  return std::get<0>(theShapes_.at(shapeType));
+}
+
+float HcalPulseShapeLookup::getTimeShift(const int shapeType) const { return std::get<1>(theShapes_.at(shapeType)); }
+
+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 std::get<2>(theShapes_.at(shapeTypes_.at(linearizedChannelNumber)));
+}
+
+const std::string& HcalPulseShapeLookup::getChannelLabel(const unsigned linearizedChannelNumber) const {
+  return std::get<0>(theShapes_.at(shapeTypes_.at(linearizedChannelNumber)));
+}
+
+float HcalPulseShapeLookup::getChannelTimeShift(const unsigned linearizedChannelNumber) const {
+  return std::get<1>(theShapes_.at(shapeTypes_.at(linearizedChannelNumber)));
+}
+
+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));
+}
+
+float HcalPulseShapeLookup::getChannelTimeShift(const DetId& id) const {
+  return getChannelTimeShift(htopo_->detId2denseId(id));
+}
+
+void HcalPulseShapeLookup::dumpToTxt(const std::string& filename, const unsigned precision) const {
+  std::ofstream of(filename.c_str());
+  if (!of.is_open()) {
+    throw cms::Exception("HcalPulseShapeLookup::dumpToTxt: file opening error")
+        << "Failed to open output file \"" << filename << '"' << std::endl;
+  }
+  if (precision)
+    of.precision(precision);
+
+  const int nShapes = theShapes_.size();
+  const unsigned nShapeTupes = shapeTypes_.size();
+  const unsigned nChannelsMapped = nShapeTupes - std::count(shapeTypes_.begin(), shapeTypes_.end(), -1);
+  of << "# nShapes = " << nShapes << " nChannelsMapped = " << nChannelsMapped << '\n';
+  for (int i = 0; i < nShapes; ++i) {
+    const std::vector<float>& data = getShape(i).data();
+    const unsigned len = data.size();
+    of << i << ' ' << getLabel(i) << ' ' << getTimeShift(i) << ' ' << len;
+    for (const float d : data)
+      of << ' ' << d;
+    of << '\n';
+  }
+
+  of << "####\n";
+  std::map<HcalSubdetector, std::string> subDetNames;
+  subDetNames[HcalBarrel] = "HB";
+  subDetNames[HcalEndcap] = "HE";
+  for (unsigned densId = 0; densId < nShapeTupes; ++densId)
+    if (shapeTypes_[densId] >= 0) {
+      const HcalDetId hcalId(htopo_->denseId2detId(densId));
+      const HcalSubdetector subdet = hcalId.subdet();
+      assert(subDetNames.find(subdet) != subDetNames.end());
+      of << subDetNames[subdet] << std::setw(4) << hcalId.ieta() << std::setw(4) << hcalId.iphi() << std::setw(3)
+         << hcalId.depth() << ' ' << shapeTypes_[densId] << '\n';
+    }
+
+  if (!of.good()) {
+    throw cms::Exception("HcalPulseShapeLookup::dumpToTxt: failed to write the pulse shapes");
+  }
+}

CalibCalorimetry/HcalPlugins/python/Hcal_Conditions_forGlobalTag_cff.py

@@ -50,6 +50,7 @@
         noiseCorrelation  = cms.vdouble(0.0),
         noiseThreshold = cms.double(0.0),
         seedThreshold  = cms.double(0.1),
+        pulseDelay = cms.double(0.0),
         doRadiationDamage = cms.bool(False)
     ),
     he = cms.PSet(
@@ -68,6 +69,7 @@
         noiseCorrelation  = cms.vdouble(0.0),
         noiseThreshold = cms.double(0.0),
         seedThreshold  = cms.double(0.1),
+        pulseDelay = cms.double(0.0),
         doRadiationDamage = cms.bool(False)
     ),
     hf = cms.PSet(
@@ -86,6 +88,7 @@
         noiseCorrelation  = cms.vdouble(0.0),
         noiseThreshold   = cms.double(0.0),
         seedThreshold    = cms.double(0.1),
+        pulseDelay = cms.double(0.0),
         doRadiationDamage = cms.bool(False)
     ),
     ho = cms.PSet(
@@ -104,6 +107,7 @@
         noiseCorrelation  = cms.vdouble(0.0),
         noiseThreshold = cms.double(0.0),
         seedThreshold  = cms.double(0.1),
+        pulseDelay = cms.double(0.0),
         doRadiationDamage = cms.bool(False)
     ),
     hbUpgrade = cms.PSet(
@@ -122,6 +126,7 @@
         noiseCorrelation  = cms.vdouble(0.26,0.254),
         noiseThreshold = cms.double(0.0),
         seedThreshold  = cms.double(0.1),
+        pulseDelay = cms.double(0.0),
         doRadiationDamage = cms.bool(True),
         radiationDamage = cms.PSet(
             temperatureBase = cms.double(20),
@@ -151,6 +156,7 @@
         noiseCorrelation  = cms.vdouble(0.26,0.254),
         noiseThreshold = cms.double(0.0),
         seedThreshold  = cms.double(0.1),
+        pulseDelay = cms.double(0.0),
         doRadiationDamage = cms.bool(True),
         radiationDamage = cms.PSet(
             temperatureBase = cms.double(20),
@@ -180,6 +186,7 @@
         noiseCorrelation  = cms.vdouble(0.0),
         noiseThreshold = cms.double(0.0),
         seedThreshold  = cms.double(0.1),
+        pulseDelay = cms.double(0.0),
         doRadiationDamage = cms.bool(False)
     ),
     # types (in order): HcalHOZecotek, HcalHOHamamatsu, HcalHEHamamatsu1, HcalHEHamamatsu2, HcalHBHamamatsu1, HcalHBHamamatsu2, HcalHPD