Adapt HGCal rechit producer for HGCal testbeam data

CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h

@@ -32,6 +32,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<>;

CondFormats/HGCalObjects/interface/HGCalMappingParameterSoA.h

@@ -61,6 +61,9 @@ namespace hgcal {
                       SOA_COLUMN(uint32_t, modInfoIdx),
                       SOA_COLUMN(uint32_t, cellInfoIdx),
                       SOA_COLUMN(uint32_t, chNumber),
+                      SOA_COLUMN(uint32_t, layer),
+                      SOA_COLUMN(float, eta),
+                      SOA_COLUMN(float, phi),
                       SOA_COLUMN(float, x),
                       SOA_COLUMN(float, y),
                       SOA_COLUMN(float, z))

Configuration/PyReleaseValidation/python/relval_standard.py

@@ -743,6 +743,8 @@
 workflows[8.1] = ['', ['BeamHalo_UP18','DIGICOS_UP18','RECOCOS_UP18','ALCABH_UP18','HARVESTCOS_UP18']]
 workflows[8.2] = ['', ['BeamHalo_UP21','DIGICOS_UP21','RECOCOS_UP21','ALCABH_UP21','HARVESTCOS_UP21']]
 
+workflows[77] = ['', ['HGCal_TestBeam']]
+
 workflows[11] = ['', ['MinBias','DIGI','RECOMIN','HARVEST','ALCAMIN']]
 workflows[28] = ['', ['QCD_Pt_80_120','DIGI','RECO','HARVEST']]
 workflows[27] = ['', ['WM','DIGI','RECO','HARVEST']]

Configuration/PyReleaseValidation/python/relval_steps.py

@@ -1292,6 +1292,11 @@ def genS(fragment,howMuch):
 # steps['CosmicsINPUT']={'INPUT':InputInfo(dataSet='/RelValCosmics/%s/GEN-SIM'%(baseDataSetRelease[0],),location='STD')}
 steps['BeamHaloINPUT']={'INPUT':InputInfo(dataSet='/RelValBeamHalo/%s/GEN-SIM'%(baseDataSetRelease[0],),location='STD')}
 
+# Phase2 HGCAL test beam
+from Configuration.PyReleaseValidation.upgradeWorkflowComponents import upgradeProperties
+hgcalTB=upgradeProperties['Run4']['Run4D104'] # so if the default changes, change wf only here
+steps['HGCal_TestBeam']={'--conditions':hgcalTB['GT'],'--era':hgcalTB['Era'],'--customise':'RecoLocalCalo/Configuration/hgcalTestBeamLocalReco_cff.runRecoForSep2024TB','-s':'NONE','--datatier':'RECO','--eventcontent':'FEVTDEBUG','--data':'','--process':'RECO','-n':'10','--filein':'/store/group/dpg_hgcal/comm_hgcal/relval/RAW2DIGI.root'}
+
 steps['QCD_Pt_50_80']=genS('QCD_Pt_50_80_8TeV_TuneCUETP8M1_cfi',Kby(25,100))
 steps['QCD_Pt_15_20']=genS('QCD_Pt_15_20_8TeV_TuneCUETP8M1_cfi',Kby(25,100))
 steps['ZTTHS']=merge([Kby(25,100),steps['ZTT']])

DataFormats/HGCalReco/interface/HGCalSoARecHits.h

@@ -14,13 +14,17 @@ GENERATE_SOA_LAYOUT(HGCalSoARecHitsLayout,
                     SOA_COLUMN(float, dim2),
                     SOA_COLUMN(float, dim3),
                     SOA_COLUMN(int, layer),
-                    SOA_COLUMN(float, weight),
+                    SOA_COLUMN(float, energy),
+                    SOA_COLUMN(float, mipEnergy),
                     SOA_COLUMN(float, sigmaNoise),
                     SOA_COLUMN(unsigned int, recHitIndex),
+                    SOA_COLUMN(uint16_t, flags),
                     SOA_COLUMN(uint32_t, detid),
                     SOA_COLUMN(float, time),
                     SOA_COLUMN(float, timeError))
 
 using HGCalSoARecHits = HGCalSoARecHitsLayout<>;
 
+enum HGCalRecHitFlags { kNormal = 0x0, kEnergyInvalid = 0x1, kTimeInvalid = 0x2 };
+
 #endif  // DataFormats_PortableTestObjects_interface_TestSoA_h

Geometry/HGCalMapping/plugins/alpaka/HGCalDenseIndexInfoESProducer.cc

@@ -114,6 +114,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
               //assign det id only for full and calibration cells
               row.detid() = 0;
+              row.layer() = 0;
               if (cell_row.t() == 1 || cell_row.t() == 0) {
                 if (isSiPM) {
                   row.detid() = ::hgcal::mappingtools::getSiPMDetId(module_row.zside(),
@@ -122,19 +123,25 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
                                                                     module_row.celltype(),
                                                                     cell_row.i1(),
                                                                     cell_row.i2());
+                  row.layer() = HGCScintillatorDetId(row.detid()).layer();
                 } else {
                   row.detid() = module_row.detid() + cell_row.detid();
+                  row.layer() = HGCSiliconDetId(row.detid()).layer();
                 }
 
                 //assign position from geometry
                 row.x() = 0;
                 row.y() = 0;
                 row.z() = 0;
+                row.eta() = 0;
+                row.phi() = 0;
                 if (hgcal_geom != nullptr) {
                   GlobalPoint position = hgcal_geom->getPosition(row.detid());
                   row.x() = position.x();
                   row.y() = position.y();
                   row.z() = position.z();
+                  row.eta() = position.eta();
+                  row.phi() = position.phi();
                 }
               }
             }  // end cell loop

Geometry/HGCalMapping/plugins/alpaka/HGCalMappingCellESProducer.cc

@@ -110,7 +110,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
             //identify special cases (Si vs SiPM, calib vs normal)
             std::string typecode = pmap.getAttr("Typecode", row);
             auto typeidx = cellIndexer.getEnumFromTypecode(typecode);
-            bool isSiPM = typecode.find("TM") != std::string::npos;
+            bool isSiPM = (typecode[0] == 'T');
             int rocpin = pmap.getIntAttr("ROCpin", row);
             int celltype = pmap.getIntAttr("t", row);
             int i1(0), i2(0), sensorcell(0);

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/Configuration/python/hgcalTestBeamLocalReco_cff.py

@@ -0,0 +1,60 @@
+import FWCore.ParameterSet.Config as cms
+
+def runRecoForSep2024TB(process):
+    process.load('Configuration.StandardSequences.Accelerators_cff')
+
+    process.load(f"Configuration.Geometry.GeometryExtendedRun4D104Reco_cff")
+    process.load(f"Configuration.Geometry.GeometryExtendedRun4D104_cff")
+    from Geometry.HGCalMapping.hgcalmapping_cff import customise_hgcalmapper
+    process = customise_hgcalmapper(process, modules='Geometry/HGCalMapping/data/ModuleMaps/modulelocator_Sep2024TBv2.txt')
+
+    # Exclude rawMetaDataCollector from input TB file
+    process.source.inputCommands = cms.untracked.vstring('keep *', 'drop *_rawMetaDataCollector_*_*')
+
+    # Keep HGCal output
+    process.RecoLocalCaloRECO.outputCommands = ['keep *_hgcal*_*_*']
+
+    # ESProducer to load calibration parameters from JSON file, like pedestals
+    process.hgcalCalibParamESProducer = cms.ESProducer('hgcalrechit::HGCalCalibrationESProducer@alpaka',
+        filename=cms.FileInPath('RecoLocalCalo/HGCalRecProducers/data/testbeam/level0_calib_Relay1727033054.json'),
+        filenameEnergyLoss=cms.FileInPath('RecoLocalCalo/HGCalRecProducers/data/testbeam/hgcal_energyloss_v16.json'),
+        indexSource=cms.ESInputTag('hgCalMappingESProducer',''),
+        mapSource=cms.ESInputTag('hgCalMappingModuleESProducer','')
+    )
+
+    process.hgcalSoARecHits = cms.EDProducer('alpaka_serial_sync::HGCalRecHitsProducer',
+        digis=cms.InputTag('hgcalDigis', ''),
+        calibSource=cms.ESInputTag('hgcalCalibParamESProducer', ''),
+        n_hits_scale=cms.int32(1),
+        n_blocks=cms.int32(1024),
+        n_threads=cms.int32(4096)
+    )
+
+    from RecoLocalCalo.HGCalRecProducers.hgCalSoARecHitsLayerClustersProducer_cfi import hgCalSoARecHitsLayerClustersProducer
+    process.hgcalSoARecHitsLayerClusters = hgCalSoARecHitsLayerClustersProducer.clone(
+        hgcalRecHitsSoA = "hgcalSoARecHits"
+    )
+
+    from RecoLocalCalo.HGCalRecProducers.hgCalSoALayerClustersProducer_cfi import hgCalSoALayerClustersProducer
+    process.hgcalSoALayerClusters = hgCalSoALayerClustersProducer.clone(
+        hgcalRecHitsLayerClustersSoA = "hgcalSoARecHitsLayerClusters",
+        hgcalRecHitsSoA = "hgcalSoARecHits"
+    )
+
+    from RecoLocalCalo.HGCalRecProducers.hgCalLayerClustersFromSoAProducer_cfi import hgCalLayerClustersFromSoAProducer
+    process.hgcalMergeLayerClusters = hgCalLayerClustersFromSoAProducer.clone(
+        hgcalRecHitsLayerClustersSoA = "hgcalSoARecHitsLayerClusters",
+        hgcalRecHitsSoA = "hgcalSoARecHits",
+        src = "hgcalSoALayerClusters"
+    )
+
+    process.reco_task = cms.Task(
+        process.hgcalSoARecHits,
+        process.hgcalSoARecHitsLayerClusters,
+        process.hgcalSoALayerClusters,
+        process.hgcalMergeLayerClusters
+    )
+    process.hgcalTestBeamLocalRecoSequence = cms.Path(process.reco_task)
+    process.schedule.insert(0, process.hgcalTestBeamLocalRecoSequence)
+
+    return process

RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h

@@ -1,6 +1,7 @@
 #ifndef RecoLocalCalo_HGCalRecAlgos_HGCalESProducerTools_h
 #define RecoLocalCalo_HGCalRecAlgos_HGCalESProducerTools_h
 
+#include "FWCore/Utilities/interface/Exception.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