Merge pull request #47829 from RSalvatico/HGCALCalibrationCell

cmsbuild · web-flow · commit b77beb911827 · 2025-05-17T19:29:12.000+02:00
Add treatment of HGCAL calibration and surrounding cells
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, caliboffset),
                       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
@@ -1,4 +1,5 @@
 #include "FWCore/ParameterSet/interface/FileInPath.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "FWCore/Utilities/interface/ESGetToken.h"
 
@@ -12,6 +13,7 @@
 #include "CondFormats/HGCalObjects/interface/HGCalMappingCellIndexer.h"
 #include "CondFormats/HGCalObjects/interface/HGCalMappingParameterHost.h"
 #include "CondFormats/HGCalObjects/interface/alpaka/HGCalMappingParameterDevice.h"
+#include "CondFormats/HGCalObjects/interface/HGCalMappingModuleIndexer.h"
 #include "DataFormats/HGCalDigi/interface/HGCalElectronicsId.h"
 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
@@ -21,6 +23,7 @@
 #include <iostream>
 #include <fstream>
 #include <sstream>
+#include <charconv>
 
 namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
@@ -30,9 +33,12 @@ 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<edm::FileInPath>("offsetfile")) {
         auto cc = setWhatProduced(this);
         cellIndexTkn_ = cc.consumes(iConfig.getParameter<edm::ESInputTag>("cellindexer"));
+        moduleIndexTkn_ = cc.consumes(iConfig.getParameter<edm::ESInputTag>("moduleindexer"));
       }
 
       //
@@ -41,18 +47,59 @@ 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<edm::ESInputTag>("moduleindexer", edm::ESInputTag(""))->setComment("Module index tool");
+        desc.add<edm::FileInPath>(
+                "offsetfile",
+                edm::FileInPath("Geometry/HGCalMapping/data/CellMaps/calibration_to_surrounding_offsetMap.txt"))
+            ->setComment("file containing the offsets between calibration and surrounding cells");
         descriptions.addWithDefaultLabel(desc);
       }
 
+      std::map<int, int> makeOffsetMap(edm::FileInPath input_offsetfile,
+                                       const HGCalMappingCellIndexer& cellIndexer,
+                                       const HGCalMappingModuleIndexer& moduleIndexer) {
+        std::map<int, int> offsetMap;
+        const auto& offsetfile = input_offsetfile.fullPath();
+        ::hgcal::mappingtools::HGCalEntityList omap;
+        edm::FileInPath fip(offsetfile);
+        omap.buildFrom(fip.fullPath());
+        auto& mapEntries = omap.getEntries();
+
+        for (auto row : mapEntries) {
+          std::string typecode = omap.getAttr("Typecode", row);
+          const auto& allTypecodes = moduleIndexer.getTypecodeMap();
+          // Skip if typecode is not in the module indexer
+          bool typecodeFound = false;
+          for (const auto& key : allTypecodes) {
+            if (key.first.find(typecode) != std::string::npos) {
+              typecodeFound = true;
+              break;
+            }
+          }
+          if (!typecodeFound)
+            continue;
+          int roc = omap.getIntAttr("ROC", row);
+          int halfroc = omap.getIntAttr("HalfROC", row);
+          int readoutsequence = omap.getIntAttr("Seq", row);
+          int offset = omap.getIntAttr("Offset", row);
+          int idx = cellIndexer.denseIndex(typecode, roc, halfroc, readoutsequence);
+          offsetMap[idx] = offset;
+        }
+        return offsetMap;
+      }
+
       //
       std::optional<HGCalMappingCellParamHost> produce(const HGCalElectronicsMappingRcd& iRecord) {
-        //get cell indexer
+        //get cell and module indexers
         const HGCalMappingCellIndexer& cellIndexer = iRecord.get(cellIndexTkn_);
+        const HGCalMappingModuleIndexer& moduleIndexer = iRecord.get(moduleIndexTkn_);
         const uint32_t size = cellIndexer.maxDenseIndex();  // channel-level size
         HGCalMappingCellParamHost cellParams(size, cms::alpakatools::host());
         for (uint32_t i = 0; i < size; i++)
           cellParams.view()[i].valid() = false;
 
+        auto offsetMap = makeOffsetMap(offsetfile_, cellIndexer, moduleIndexer);
+
         //loop over cell types and then over cells
         for (const auto& url : filelist_) {
           ::hgcal::mappingtools::HGCalEntityList pmap;
@@ -113,6 +160,10 @@ 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;
+
+            int offset = (iscalib && offsetMap.find(idx) != offsetMap.end()) ? offsetMap[idx] : 0;
+            cell.caliboffset() = offset;
+
           }  //end loop over entities
         }  //end loop over cell types
 
@@ -121,7 +172,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
     private:
       edm::ESGetToken<HGCalMappingCellIndexer, HGCalElectronicsMappingRcd> cellIndexTkn_;
+      edm::ESGetToken<HGCalMappingModuleIndexer, HGCalElectronicsMappingRcd> moduleIndexTkn_;
       const std::vector<std::string> filelist_;
+      edm::FileInPath 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/calibration_to_surrounding_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.FileInPath(offsetfile))
     process.hgCalMappingModuleESProducer = cms.ESProducer('hgcal::HGCalMappingModuleESProducer@alpaka',
                                                           filename=cms.FileInPath(modules),
                                                           moduleindexer=cms.ESInputTag(''))
diff --git a/Geometry/HGCalMapping/test/testMappingModuleIndexer_cfg.py b/Geometry/HGCalMapping/test/testMappingModuleIndexer_cfg.py
@@ -9,6 +9,9 @@
                  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,
+                 info="Path to calibration-to-surrounding cell offset file. Absolute, or relative to CMSSW src directory")
+
 options.parseArguments()
 
 process.source = cms.Source('EmptySource')
@@ -22,7 +25,8 @@
 process = customise_hgcalmapper(process,
                                 modules=options.modules,
                                 sicells=options.sicells,
-                                sipmcells=options.sipmcells)
+                                sipmcells=options.sipmcells,
+                                offsetfile=options.offsetfile)
 
 #Geometry
 process.load('Configuration.Geometry.GeometryExtended2026D99Reco_cff')
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].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
+        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
