add treatment of HGCAL calibration and surrounding cells

RSalvatico · RSalvatico · commit 47143af61565 · 2025-04-09T15:48:37.000+02:00
diff --git a/CondFormats/HGCalObjects/interface/HGCalMappingParameterSoA.h b/CondFormats/HGCalObjects/interface/HGCalMappingParameterSoA.h
@@ -35,6 +35,7 @@ namespace hgcal {
                       SOA_COLUMN(bool, valid),
                       SOA_COLUMN(bool, isHD),
                       SOA_COLUMN(bool, iscalib),
+                      SOA_COLUMN(int, offset),
                       SOA_COLUMN(bool, isSiPM),
                       SOA_COLUMN(uint16_t, typeidx),
                       SOA_COLUMN(uint16_t, chip),
diff --git a/Geometry/HGCalMapping/plugins/alpaka/HGCalMappingCellESProducer.cc b/Geometry/HGCalMapping/plugins/alpaka/HGCalMappingCellESProducer.cc
@@ -30,7 +30,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     public:
       //
       HGCalMappingCellESProducer(const edm::ParameterSet& iConfig)
-          : ESProducer(iConfig), filelist_(iConfig.getParameter<std::vector<std::string> >("filelist")) {
+          : ESProducer(iConfig),
+            filelist_(iConfig.getParameter<std::vector<std::string> >("filelist")),
+            offsetfile_(iConfig.getParameter<std::string>("offsetfile")) {
         auto cc = setWhatProduced(this);
         cellIndexTkn_ = cc.consumes(iConfig.getParameter<edm::ESInputTag>("cellindexer"));
       }
@@ -41,9 +43,41 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         desc.add<std::vector<std::string> >("filelist", std::vector<std::string>({}))
             ->setComment("list of files with the readout cells of each module");
         desc.add<edm::ESInputTag>("cellindexer", edm::ESInputTag(""))->setComment("Dense cell index tool");
+        desc.add<std::string>("offsetfile", std::string({}))
+            ->setComment("file containing the offsets between calibration and surrounding cells");
         descriptions.addWithDefaultLabel(desc);
       }
 
+      std::map<std::tuple<std::string, int, int>, int> makeOffsetMap(std::string input_offsetfile) {
+        std::map<std::tuple<std::string, int, int>, int> offsetMap;
+        auto offsetfile = edm::FileInPath(input_offsetfile).fullPath();
+        std::ifstream stream_offsetfile(offsetfile);
+        std::string line;
+        // Skip the first line (description of column content)
+        std::getline(stream_offsetfile, line);
+        while (std::getline(stream_offsetfile, line)) {
+          std::istringstream iss(line);
+          std::string r_typecode_str, r_chip_str, r_half_str, offset_str;
+          if (!(iss >> r_typecode_str >> r_chip_str >> r_half_str >> offset_str)) {
+            std::cout << "Error reading offset file" << std::endl;
+            break;
+          }
+          try {
+            int r_chip = std::stoi(r_chip_str);
+            int r_half = std::stoi(r_half_str);
+            int offset = std::stoi(offset_str);
+            offsetMap[std::make_tuple(r_typecode_str, r_chip, r_half)] = offset;
+          } catch (const std::invalid_argument& e) {
+            std::cout << "Invalid argument: " << e.what() << std::endl;
+            break;
+          } catch (const std::out_of_range& e) {
+            std::cout << "Out of range: " << e.what() << std::endl;
+            break;
+          }
+        }
+        return offsetMap;
+      }
+
       //
       std::optional<HGCalMappingCellParamHost> produce(const HGCalElectronicsMappingRcd& iRecord) {
         //get cell indexer
@@ -53,6 +87,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         for (uint32_t i = 0; i < size; i++)
           cellParams.view()[i].valid() = false;
 
+        auto offsetMap = makeOffsetMap(offsetfile_);
+
         //loop over cell types and then over cells
         for (const auto& url : filelist_) {
           ::hgcal::mappingtools::HGCalEntityList pmap;
@@ -113,6 +149,12 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
             cell.trace() = pmap.getFloatAttr("trace", row);
             cell.eleid() = HGCalElectronicsId(false, 0, 0, 0, chip * 2 + half, seq).raw();
             cell.detid() = detid;
+
+            //calibration cell-to-surrounding cell offset. Can be !=0 only for calibration cells
+            auto mapKey = std::make_tuple(typecode, chip, half);
+            int offset = (iscalib && offsetMap.find(mapKey) != offsetMap.end()) ? offsetMap[mapKey] : 0;
+            cell.offset() = offset;
+
           }  //end loop over entities
         }  //end loop over cell types
 
@@ -122,6 +164,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     private:
       edm::ESGetToken<HGCalMappingCellIndexer, HGCalElectronicsMappingRcd> cellIndexTkn_;
       const std::vector<std::string> filelist_;
+      std::string offsetfile_;
     };
 
   }  // namespace hgcal
diff --git a/Geometry/HGCalMapping/python/hgcalmapping_cff.py b/Geometry/HGCalMapping/python/hgcalmapping_cff.py
@@ -3,7 +3,8 @@
 def customise_hgcalmapper(process,
                           modules = 'Geometry/HGCalMapping/data/ModuleMaps/modulelocator_test.txt',
                           sicells = 'Geometry/HGCalMapping/data/CellMaps/WaferCellMapTraces.txt',
-                          sipmcells = 'Geometry/HGCalMapping/data/CellMaps/channels_sipmontile.hgcal.txt'):
+                          sipmcells = 'Geometry/HGCalMapping/data/CellMaps/channels_sipmontile.hgcal.txt',
+                          offsetfile = 'Geometry/HGCalMapping/data/CellMaps/offsetMap.txt'):
     """the following function configures the mapping producers
     NOTE: for production-targetted configs should be avoided as it checks if the process as 
     already the Accelerators sequence loaded, if not it loads it to the process"""
@@ -18,7 +19,8 @@ def customise_hgcalmapper(process,
         
     process.hgCalMappingCellESProducer = cms.ESProducer('hgcal::HGCalMappingCellESProducer@alpaka',
                                                         filelist=cms.vstring(sicells, sipmcells),
-                                                        cellindexer=cms.ESInputTag(''))
+                                                        cellindexer=cms.ESInputTag(''),
+                                                        offsetfile=cms.string(offsetfile))
     process.hgCalMappingModuleESProducer = cms.ESProducer('hgcal::HGCalMappingModuleESProducer@alpaka',
                                                           filename=cms.FileInPath(modules),
                                                           moduleindexer=cms.ESInputTag(''))
diff --git a/RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h b/RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h
@@ -16,11 +16,14 @@
 #include "DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h"
 #include "CondFormats/HGCalObjects/interface/HGCalCalibrationParameterHost.h"
 #include "CondFormats/HGCalObjects/interface/alpaka/HGCalCalibrationParameterDevice.h"
+#include "CondFormats/HGCalObjects/interface/HGCalMappingParameterHost.h"
+#include "CondFormats/HGCalObjects/interface/alpaka/HGCalMappingParameterDevice.h"
 
 namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
   using namespace hgcaldigi;
   using namespace hgcalrechit;
+  using namespace hgcal;
 
   class HGCalRecHitCalibrationAlgorithms {
   public:
@@ -29,7 +32,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     HGCalRecHitDevice calibrate(Queue& queue,
                                 HGCalDigiHost const& host_digis,
                                 HGCalCalibParamDevice const& device_calib,
-                                HGCalConfigParamDevice const& device_config) const;
+                                HGCalConfigParamDevice const& device_config,
+                                HGCalMappingCellParamDevice const& device_mapping,
+                                HGCalDenseIndexInfoDevice const& device_index) const;
 
   private:
     void print(HGCalDigiHost const& digis, int max = -1) const;
diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc
@@ -28,9 +28,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         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) && (digiflags != hgcal::DIGI_FLAG::ZS_ToA_ADCm1));
+        bool isAvailable((digiflags != ::hgcal::DIGI_FLAG::Invalid) &&
+                         (digiflags != ::hgcal::DIGI_FLAG::NotAvailable) && calibvalid);
+        bool isToAavailable((digiflags != ::hgcal::DIGI_FLAG::ZS_ToA) &&
+                            (digiflags != ::hgcal::DIGI_FLAG::ZS_ToA_ADCm1));
         recHits[idx].flags() = (!isAvailable) * hgcalrechit::HGCalRecHitFlags::EnergyInvalid +
                                (!isToAavailable) * hgcalrechit::HGCalRecHitFlags::TimeInvalid;
       }
@@ -61,8 +62,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         bool calibvalid = calib.valid();
         auto digi = digis[idx];
         auto digiflags = digi.flags();
-        bool isAvailable((digiflags != hgcal::DIGI_FLAG::Invalid) && (digiflags != hgcal::DIGI_FLAG::NotAvailable) &&
-                         calibvalid);
+        bool isAvailable((digiflags != ::hgcal::DIGI_FLAG::Invalid) &&
+                         (digiflags != ::hgcal::DIGI_FLAG::NotAvailable) && calibvalid);
         bool useTOT((digi.tctp() == 3) && isAvailable);
         bool useADC(!useTOT && isAvailable);
         recHits[idx].energy() = useADC * adc_denoise(digi.adc(),
@@ -111,9 +112,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         bool calibvalid = calib.valid();
         auto digi = digis[idx];
         auto digiflags = digi.flags();
-        bool isAvailable((digiflags != hgcal::DIGI_FLAG::Invalid) && (digiflags != hgcal::DIGI_FLAG::NotAvailable) &&
-                         calibvalid);
-        bool isToAavailable((digiflags != hgcal::DIGI_FLAG::ZS_ToA) && (digiflags != hgcal::DIGI_FLAG::ZS_ToA_ADCm1));
+        bool isAvailable((digiflags != ::hgcal::DIGI_FLAG::Invalid) &&
+                         (digiflags != ::hgcal::DIGI_FLAG::NotAvailable) && calibvalid);
+        bool isToAavailable((digiflags != ::hgcal::DIGI_FLAG::ZS_ToA) &&
+                            (digiflags != ::hgcal::DIGI_FLAG::ZS_ToA_ADCm1));
         bool isGood(isAvailable && isToAavailable);
         //INL correction
         auto toa = isGood * toa_inl_corr(digi.toa(), calib.TOA_CTDC(), calib.TOA_FTDC());
@@ -125,6 +127,59 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     }
   };
 
+  struct HGCalRecHitCalibrationKernel_handleCalibCell {
+    ALPAKA_FN_ACC void operator()(Acc1D const& acc,
+                                  HGCalDigiDevice::View digis,
+                                  HGCalRecHitDevice::View recHits,
+                                  HGCalCalibParamDevice::ConstView calibs,
+                                  HGCalMappingCellParamDevice::ConstView maps,
+                                  HGCalDenseIndexInfoDevice::ConstView index) const {
+      auto time_average = [&](float time_surr, float time_calib, float energy_surr, float energy_calib) {
+        bool is_time_surr(time_surr > 0);
+        bool is_time_calib(time_calib > 0);
+        float totalEn = (is_time_surr * energy_surr + is_time_calib * energy_calib);
+        float weighted_average =
+            (totalEn > 0)
+                ? (is_time_surr * energy_surr * time_surr + is_time_calib * energy_calib * time_calib) / totalEn
+                : 0.0f;
+        return weighted_average;
+      };
+
+      for (auto idx : uniform_elements(acc, digis.metadata().size())) {
+        auto calib = calibs[idx];
+        bool calibvalid = calib.valid();
+        auto digi = digis[idx];
+        auto digiflags = digi.flags();
+        bool isAvailable((digiflags != ::hgcal::DIGI_FLAG::Invalid) &&
+                         (digiflags != ::hgcal::DIGI_FLAG::NotAvailable) && calibvalid);
+        bool isToAavailable((digiflags != ::hgcal::DIGI_FLAG::ZS_ToA) &&
+                            (digiflags != ::hgcal::DIGI_FLAG::ZS_ToA_ADCm1));
+
+        auto cellIndex = index[idx].cellInfoIdx();
+        bool isCalibCell(maps[cellIndex].iscalib());
+        int offset = maps[cellIndex].offset();
+        bool is_surr_cell((offset != 0) && isAvailable && isCalibCell);
+
+        //Effectively operate only on the cell that surrounds the calibration cells
+        if (!is_surr_cell) {
+          continue;
+        }
+
+        recHits[idx + offset].flags() = hgcalrechit::HGCalRecHitFlags::Normal;
+
+        recHits[idx + offset].time() = isToAavailable * time_average(recHits[idx + offset].time(),
+                                                                     recHits[idx].time(),
+                                                                     recHits[idx + offset].energy(),
+                                                                     recHits[idx].energy());
+
+        bool is_negative_surr_energy(recHits[idx + offset].energy() < 0);
+        auto negative_energy_correction = (-1.0 * recHits[idx + offset].energy()) * is_negative_surr_energy;
+
+        recHits[idx + offset].energy() += (negative_energy_correction + recHits[idx].energy());
+      }
+    }
+  };
+
   struct HGCalRecHitCalibrationKernel_printRecHits {
     ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalRecHitDevice::ConstView view, int size) const {
       for (int i = 0; i < size; ++i) {
@@ -137,7 +192,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
   HGCalRecHitDevice HGCalRecHitCalibrationAlgorithms::calibrate(Queue& queue,
                                                                 HGCalDigiHost const& host_digis,
                                                                 HGCalCalibParamDevice const& device_calib,
-                                                                HGCalConfigParamDevice const& device_config) const {
+                                                                HGCalConfigParamDevice const& device_config,
+                                                                HGCalMappingCellParamDevice const& device_mapping,
+                                                                HGCalDenseIndexInfoDevice const& device_index) const {
     LogDebug("HGCalRecHitCalibrationAlgorithms") << "\n\nINFO -- Start of calibrate\n\n" << std::endl;
 
     LogDebug("HGCalRecHitCalibrationAlgorithms") << "\n\nINFO -- Copying the digis to the device\n\n" << std::endl;
@@ -176,6 +233,14 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
                         device_digis.view(),
                         device_recHits.view(),
                         device_calib.view());
+    alpaka::exec<Acc1D>(queue,
+                        grid,
+                        HGCalRecHitCalibrationKernel_handleCalibCell{},
+                        device_digis.view(),
+                        device_recHits.view(),
+                        device_calib.view(),
+                        device_mapping.view(),
+                        device_index.view());
 
     LogDebug("HGCalRecHitCalibrationAlgorithms") << "Input recHits: " << std::endl;
 #ifdef EDM_ML_DEBUG
diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitProducers.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitProducers.cc
@@ -32,6 +32,8 @@
 #include "CondFormats/HGCalObjects/interface/HGCalCalibrationParameterHost.h"
 #include "CondFormats/HGCalObjects/interface/alpaka/HGCalCalibrationParameterDevice.h"
 #include "RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h"
+#include "CondFormats/HGCalObjects/interface/alpaka/HGCalMappingParameterDevice.h"
+#include "CondFormats/DataRecord/interface/HGCalDenseIndexInfoRcd.h"
 
 // flag to assist the computational performance test
 // #define HGCAL_PERF_TEST
@@ -52,6 +54,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     const edm::EDGetTokenT<hgcaldigi::HGCalDigiHost> digisToken_;
     edm::ESGetToken<hgcalrechit::HGCalCalibParamHost, HGCalModuleConfigurationRcd> calibToken_;
     device::ESGetToken<hgcalrechit::HGCalConfigParamDevice, HGCalModuleConfigurationRcd> configToken_;
+    device::ESGetToken<hgcal::HGCalMappingCellParamDevice, HGCalElectronicsMappingRcd> mappingToken_;
+    device::ESGetToken<hgcal::HGCalDenseIndexInfoDevice, HGCalDenseIndexInfoRcd> indexingToken_;
     const device::EDPutToken<hgcalrechit::HGCalRecHitDevice> recHitsToken_;
     const HGCalRecHitCalibrationAlgorithms calibrator_;
     const int n_hits_scale;
@@ -62,6 +66,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
         digisToken_{consumes<hgcaldigi::HGCalDigiHost>(iConfig.getParameter<edm::InputTag>("digis"))},
         calibToken_{esConsumes(iConfig.getParameter<edm::ESInputTag>("calibSource"))},
         configToken_{esConsumes(iConfig.getParameter<edm::ESInputTag>("configSource"))},
+        mappingToken_{esConsumes(iConfig.getParameter<edm::ESInputTag>("mappingSource"))},
+        indexingToken_{esConsumes(iConfig.getParameter<edm::ESInputTag>("indexingSource"))},
         recHitsToken_{produces()},
         calibrator_{iConfig.getParameter<int>("n_blocks"), iConfig.getParameter<int>("n_threads")},
         n_hits_scale{iConfig.getParameter<int>("n_hits_scale")} {
@@ -77,6 +83,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     desc.add<edm::InputTag>("digis", edm::InputTag("hgcalDigis", "DIGI", "TEST"));
     desc.add("calibSource", edm::ESInputTag{})->setComment("Label for calibration parameters");
     desc.add("configSource", edm::ESInputTag{})->setComment("Label for ROC configuration parameters");
+    desc.add("mappingSource", edm::ESInputTag{})->setComment("Label for cell mapping parameters");
+    desc.add("indexingSource", edm::ESInputTag{})->setComment("Label for cell dense indexer");
     desc.add<int>("n_blocks", -1);
     desc.add<int>("n_threads", -1);
     desc.add<int>("n_hits_scale", -1);
@@ -89,6 +97,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     // Read digis
     auto const& hostCalibParamProvider = iSetup.getData(calibToken_);
     auto const& deviceConfigParamProvider = iSetup.getData(configToken_);
+    auto const& deviceMappingCellParamProvider = iSetup.getData(mappingToken_);
+    auto const& deviceIndexingParamProvider = iSetup.getData(indexingToken_);
     auto const& hostDigisIn = iEvent.get(digisToken_);
 
     //printout new conditions if available
@@ -161,7 +171,12 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     HGCalRecHitDevice recHits(oldSize, queue);
     alpaka::memcpy(queue, recHits.buffer(), tmpRecHits.const_buffer(), oldSize);
 #else
-    auto recHits = calibrator_.calibrate(queue, hostDigis, deviceCalibParam, deviceConfigParamProvider);
+    auto recHits = calibrator_.calibrate(queue,
+                                         hostDigis,
+                                         deviceCalibParam,
+                                         deviceConfigParamProvider,
+                                         deviceMappingCellParamProvider,
+                                         deviceIndexingParamProvider);
 #endif
 
 #ifdef EDM_ML_DEBUG
