[HGCAL] Add energy loss to Calibration SoA

CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h

@@ -33,6 +33,7 @@ namespace hgcalrechit {
                       SOA_EIGEN_COLUMN(Vector8f, TOA_FTDC),   // TOA fine TDC correction
                       SOA_EIGEN_COLUMN(Vector3f, TOA_TW),     // TOA timewalk correction
                       SOA_COLUMN(float, MIPS_scale),          // MIPS scale
+                      SOA_COLUMN(float, EM_scale),            // electromagnetic scale
                       SOA_COLUMN(unsigned char, valid)        // only 1 bit used: if false = mask dead channel
   )
   using HGCalCalibParamSoA = HGCalCalibParamSoALayout<>;

Geometry/HGCalMapping/test/testMappingModuleIndexer_cfg.py

@@ -9,7 +9,7 @@
                  info="Path to Si cell mapper. Absolute, or relative to CMSSW src directory")
 options.register('sipmcells','Geometry/HGCalMapping/data/CellMaps/channels_sipmontile.hgcal.txt',mytype=VarParsing.varType.string,
                  info="Path to SiPM-on-tile cell mapper. Absolute, or relative to CMSSW src directory")
-options.register('offsetfile','Geometry/HGCalMapping/data/CellMaps/calibration_to_surrounding_offsetMap.txt',mytype=VarParsing.varType.cms.FileInPath,
+options.register('offsetfile','Geometry/HGCalMapping/data/CellMaps/calibration_to_surrounding_offsetMap.txt',mytype=VarParsing.varType.string,
                  info="Path to calibration-to-surrounding cell offset file. Absolute, or relative to CMSSW src directory")
 
 options.parseArguments()
@@ -20,18 +20,15 @@
     input = cms.untracked.int32(1)
 )
 
-#electronics mapping
+# electronics mapping
 from Geometry.HGCalMapping.hgcalmapping_cff import customise_hgcalmapper
-process = customise_hgcalmapper(process,
-                                modules=options.modules,
-                                sicells=options.sicells,
-                                sipmcells=options.sipmcells,
-                                offsetfile=options.offsetfile)
+kwargs = { k: getattr(options,k) for k in ['modules','sicells','sipmcells','offsetfile'] if getattr(options,k)!='' }
+process = customise_hgcalmapper(process, **kwargs)
 
-#Geometry
-process.load('Configuration.Geometry.GeometryExtended2026D99Reco_cff')
+# Geometry
+process.load('Configuration.Geometry.GeometryExtendedRun4D104Reco_cff')
 
-#tester
+# tester
 process.tester = cms.EDAnalyzer('HGCalMappingESSourceTester')
 
 process.p = cms.Path(process.tester)

RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h

@@ -1,6 +1,7 @@
 #ifndef RecoLocalCalo_HGCalRecAlgos_HGCalESProducerTools_h
 #define RecoLocalCalo_HGCalRecAlgos_HGCalESProducerTools_h
 
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include <string>
 #include <nlohmann/json.hpp>
 using json = nlohmann::ordered_json;  // ordered_json preserves key insertion order
@@ -13,6 +14,57 @@ namespace hgcal {
                   const std::string& firstkey,
                   const std::vector<std::string>& keys,
                   const std::string& fname);
+  bool check_keys(const json& data, const std::vector<std::string>& keys, const std::string& fname);
+
+  // @short fill full SoA column with single value
+  template <typename T>
+  void fill_SoA_column_single(T* column_SoA, const float& value, const int offset, const int nrows) {
+    std::fill(column_SoA + offset, column_SoA + offset + nrows, value);
+  }
+
+  // @short fill SoA column with data from vector for any type with some offset
+  template <typename T>
+  void fill_SoA_column(
+      T* column_SoA, const std::vector<T>& values, const int offset, const int nrows, int arr_offset = 0) {
+    const int nrows_vals = values.size();
+    if (arr_offset < 0) {
+      arr_offset = 0;
+      if (nrows_vals != arr_offset + nrows) {
+        cms::Exception ex("InvalidData");
+        ex << " Expected " << nrows << " rows, but got " << nrows_vals << "!";
+        ex.addContext("Calling hgcal::fill_SoA_eigen_row()");
+        throw ex;
+      }
+    } else if (nrows_vals < arr_offset + nrows) {
+      cms::Exception ex("InvalidData");
+      ex << " Tried to copy " << nrows << " rows to SoA with offset " << arr_offset << ", but only have " << nrows_vals
+         << " values in JSON!";
+      ex.addContext("Calling hgcal::fill_SoA_eigen_row()");
+      throw ex;
+    }
+    auto begin = values.begin() + arr_offset;
+    auto end = (begin + nrows > values.end()) ? values.end() : begin + nrows;
+    std::copy(begin, end, &column_SoA[offset]);
+  }
+
+  // @short fill full SoA column with data from vector for any type
+  template <typename T, typename P>
+  void fill_SoA_eigen_row(P& soa, const std::vector<std::vector<T>>& values, const size_t row) {
+    if (row >= values.size()) {
+      cms::Exception ex("InvalidData");
+      ex << " Tried to copy row " << row << " to SoA, but only have " << values.size() << " values in JSON!";
+      ex.addContext("Calling hgcal::fill_SoA_eigen_row()");
+      throw ex;
+    }
+    if (!values.empty() && int(values[row].size()) != soa.size()) {
+      cms::Exception ex("InvalidData");
+      ex << " Expected " << soa.size() << " elements in Eigen vector, but got " << values[row].size() << "!";
+      ex.addContext("Calling hgcal::fill_SoA_eigen_row()");
+      throw ex;
+    }
+    for (int i = 0; i < soa.size(); i++)
+      soa(i) = values[row][i];
+  }
 
 }  // namespace hgcal
 

RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc

@@ -16,7 +16,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
   using namespace cms::alpakatools;
 
-  //
+  // @short set RecHit flags based on DIGI flags (data corruption, zero-suppression)
   struct HGCalRecHitCalibrationKernel_flagRecHits {
     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
                                   HGCalDigiDevice::View digis,
@@ -27,7 +27,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         bool calibvalid = calib.valid();
         auto digi = digis[idx];
         auto digiflags = digi.flags();
-        //recHits[idx].flags() = digiflags;
         bool isAvailable((digiflags != ::hgcal::DIGI_FLAG::Invalid) &&
                          (digiflags != ::hgcal::DIGI_FLAG::NotAvailable) && calibvalid);
         bool isToAavailable((digiflags != ::hgcal::DIGI_FLAG::ZS_ToA) &&
@@ -38,8 +37,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     }
   };
 
-  //
-  struct HGCalRecHitCalibrationKernel_adcToCharge {
+  // @short subtract pedestals, linearize ADC & TOT to charge (fC), and convert to energy (GeV)
+  struct HGCalRecHitCalibrationKernel_adcToEnergy {
     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
                                   HGCalDigiDevice::View digis,
                                   HGCalRecHitDevice::View recHits,
@@ -81,13 +80,13 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
                                                            calib.TOT_P1(),
                                                            calib.TOT_P2());
 
-        //after denoising/linearization apply the MIP scale
-        recHits[idx].energy() *= calib.MIPS_scale();
+        // after denoising/linearization apply the MIP & EM scale to convert to energy (GeV)
+        recHits[idx].energy() *= calib.MIPS_scale() * calib.EM_scale();
       }
     }
   };
 
-  //
+  // @short apply INL and timewalk corrections, and convert TOA to time (ps)
   struct HGCalRecHitCalibrationKernel_toaToTime {
     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
                                   HGCalDigiDevice::View digis,
@@ -117,16 +116,18 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         bool isToAavailable((digiflags != ::hgcal::DIGI_FLAG::ZS_ToA) &&
                             (digiflags != ::hgcal::DIGI_FLAG::ZS_ToA_ADCm1));
         bool isGood(isAvailable && isToAavailable);
-        //INL correction
+        // INL correction
         auto toa = isGood * toa_inl_corr(digi.toa(), calib.TOA_CTDC(), calib.TOA_FTDC());
-        //timewalk correction
+        // timewalk correction
         toa = isGood * toa_tw_corr(toa, recHits[idx].energy(), calib.TOA_TW());
-        //toa to ps
+        // toa to ps
         recHits[idx].time() = toa * hgcalrechit::TOAtops;
       }
     }
   };
 
+  // @short sum the energy of the RecHits in the calibration and surrounding cells
+  // see https://github.com/cms-sw/cmssw/pull/47829
   struct HGCalRecHitCalibrationKernel_handleCalibCell {
     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
                                   HGCalDigiDevice::View digis,
@@ -157,10 +158,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
         auto cellIndex = index[idx].cellInfoIdx();
         bool isCalibCell(maps[cellIndex].iscalib());
-        int offset = maps[cellIndex].caliboffset();  //Calibration-to-surrounding cell offset
+        int offset = maps[cellIndex].caliboffset();  // calibration-to-surrounding cell offset
         bool is_surr_cell((offset != 0) && isAvailable && isCalibCell);
 
-        //Effectively operate only on the cell that surrounds the calibration cells
+        // effectively operate only on the cell that surrounds the calibration cells
         if (!is_surr_cell) {
           continue;
         }
@@ -180,6 +181,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     }
   };
 
+  // @short print RecHits for debugging
   struct HGCalRecHitCalibrationKernel_printRecHits {
     ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalRecHitDevice::ConstView view, int size) const {
       for (int i = 0; i < size; ++i) {
@@ -189,6 +191,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     }
   };
 
+  // @short apply all calibration kernels
   HGCalRecHitDevice HGCalRecHitCalibrationAlgorithms::calibrate(Queue& queue,
                                                                 HGCalDigiHost const& host_digis,
                                                                 HGCalCalibParamDevice const& device_calib,
@@ -222,7 +225,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
                         device_calib.view());
     alpaka::exec<Acc1D>(queue,
                         grid,
-                        HGCalRecHitCalibrationKernel_adcToCharge{},
+                        HGCalRecHitCalibrationKernel_adcToEnergy{},
                         device_digis.view(),
                         device_recHits.view(),
                         device_calib.view());