diff --git a/CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h b/CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h index 63c8aa844a748..fce76b38bc722 100644 --- a/CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h +++ b/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<>; diff --git a/CondFormats/HGCalObjects/interface/HGCalMappingParameterSoA.h b/CondFormats/HGCalObjects/interface/HGCalMappingParameterSoA.h index ef7e6500b5a35..1574d778fd1c3 100644 --- a/CondFormats/HGCalObjects/interface/HGCalMappingParameterSoA.h +++ b/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)) diff --git a/Configuration/PyReleaseValidation/python/relval_standard.py b/Configuration/PyReleaseValidation/python/relval_standard.py index 626f49c04a48f..566a85a51ed4f 100644 --- a/Configuration/PyReleaseValidation/python/relval_standard.py +++ b/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']] diff --git a/Configuration/PyReleaseValidation/python/relval_steps.py b/Configuration/PyReleaseValidation/python/relval_steps.py index 1c757c038437d..6d95c6d6d10ae 100644 --- a/Configuration/PyReleaseValidation/python/relval_steps.py +++ b/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']]) diff --git a/DataFormats/HGCalRecHit/BuildFile.xml b/DataFormats/HGCalRecHit/BuildFile.xml deleted file mode 100644 index 543d2de5864fc..0000000000000 --- a/DataFormats/HGCalRecHit/BuildFile.xml +++ /dev/null @@ -1,10 +0,0 @@ - - - - - - - - - - diff --git a/DataFormats/HGCalRecHit/interface/HGCalRecHitHost.h b/DataFormats/HGCalRecHit/interface/HGCalRecHitHost.h deleted file mode 100644 index 1a29f8f54c8d1..0000000000000 --- a/DataFormats/HGCalRecHit/interface/HGCalRecHitHost.h +++ /dev/null @@ -1,14 +0,0 @@ -#ifndef DataFormats_HGCalRecHit_interface_HGCalRecHitHost_h -#define DataFormats_HGCalRecHit_interface_HGCalRecHitHost_h - -#include "DataFormats/Portable/interface/PortableHostCollection.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h" - -namespace hgcalrechit { - - // SoA with x, y, z, id fields in host memory - using HGCalRecHitHost = PortableHostCollection; - -} // namespace hgcalrechit - -#endif // DataFormats_HGCalRecHit_interface_HGCalRecHitHost_h diff --git a/DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h b/DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h deleted file mode 100644 index 87053091a88a1..0000000000000 --- a/DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h +++ /dev/null @@ -1,20 +0,0 @@ -#ifndef DataFormats_HGCalRecHit_interface_HGCalRecHitSoA_h -#define DataFormats_HGCalRecHit_interface_HGCalRecHitSoA_h - -#include "DataFormats/SoATemplate/interface/SoACommon.h" -#include "DataFormats/SoATemplate/interface/SoALayout.h" - -namespace hgcalrechit { - - // Generate structure of arrays (SoA) layout with RecHit dataformat - GENERATE_SOA_LAYOUT(HGCalRecHitSoALayout, - SOA_COLUMN(double, energy), - SOA_COLUMN(double, time), - SOA_COLUMN(uint16_t, flags)) - using HGCalRecHitSoA = HGCalRecHitSoALayout<>; - - enum HGCalRecHitFlags { Normal = 0x0, EnergyInvalid = 0x1, TimeInvalid = 0x2 }; - -} // namespace hgcalrechit - -#endif // DataFormats_HGCalRecHit_interface_HGCalRecHitSoA_h diff --git a/DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h b/DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h deleted file mode 100644 index 158b8a1f90107..0000000000000 --- a/DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h +++ /dev/null @@ -1,23 +0,0 @@ -#ifndef DataFormats_HGCalRecHit_interface_alpaka_HGCalRecHitDevice_h -#define DataFormats_HGCalRecHit_interface_alpaka_HGCalRecHitDevice_h - -#include "DataFormats/Portable/interface/alpaka/PortableCollection.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h" -#include "HeterogeneousCore/AlpakaInterface/interface/config.h" - -namespace ALPAKA_ACCELERATOR_NAMESPACE { - - namespace hgcalrechit { - - // make the names from the top-level hgcalrechit namespace visible for unqualified lookup - // inside the ALPAKA_ACCELERATOR_NAMESPACE::hgcalrechit namespace - using namespace ::hgcalrechit; - - // SoA in device global memory - using HGCalRecHitDevice = PortableCollection; - - } // namespace hgcalrechit - -} // namespace ALPAKA_ACCELERATOR_NAMESPACE - -#endif // DataFormats_HGCalRecHit_interface_alpaka_HGCalRecHitDevice_h diff --git a/DataFormats/HGCalRecHit/src/alpaka/classes_cuda.h b/DataFormats/HGCalRecHit/src/alpaka/classes_cuda.h deleted file mode 100644 index 0f81a12c9ec38..0000000000000 --- a/DataFormats/HGCalRecHit/src/alpaka/classes_cuda.h +++ /dev/null @@ -1,4 +0,0 @@ -#include "DataFormats/Common/interface/DeviceProduct.h" -#include "DataFormats/Common/interface/Wrapper.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h" -#include "DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h" diff --git a/DataFormats/HGCalRecHit/src/alpaka/classes_cuda_def.xml b/DataFormats/HGCalRecHit/src/alpaka/classes_cuda_def.xml deleted file mode 100644 index aa9fef56baae9..0000000000000 --- a/DataFormats/HGCalRecHit/src/alpaka/classes_cuda_def.xml +++ /dev/null @@ -1,5 +0,0 @@ - - - - - diff --git a/DataFormats/HGCalRecHit/src/alpaka/classes_rocm.h b/DataFormats/HGCalRecHit/src/alpaka/classes_rocm.h deleted file mode 100644 index 0f81a12c9ec38..0000000000000 --- a/DataFormats/HGCalRecHit/src/alpaka/classes_rocm.h +++ /dev/null @@ -1,4 +0,0 @@ -#include "DataFormats/Common/interface/DeviceProduct.h" -#include "DataFormats/Common/interface/Wrapper.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h" -#include "DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h" diff --git a/DataFormats/HGCalRecHit/src/alpaka/classes_rocm_def.xml b/DataFormats/HGCalRecHit/src/alpaka/classes_rocm_def.xml deleted file mode 100644 index e3bb4ace8edf0..0000000000000 --- a/DataFormats/HGCalRecHit/src/alpaka/classes_rocm_def.xml +++ /dev/null @@ -1,5 +0,0 @@ - - - - - diff --git a/DataFormats/HGCalRecHit/src/classes.cc b/DataFormats/HGCalRecHit/src/classes.cc deleted file mode 100644 index b7e3917ad8ceb..0000000000000 --- a/DataFormats/HGCalRecHit/src/classes.cc +++ /dev/null @@ -1,4 +0,0 @@ -#include "DataFormats/Portable/interface/PortableHostCollectionReadRules.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitHost.h" - -SET_PORTABLEHOSTCOLLECTION_READ_RULES(hgcalrechit::HGCalRecHitHost); diff --git a/DataFormats/HGCalRecHit/src/classes.h b/DataFormats/HGCalRecHit/src/classes.h deleted file mode 100644 index 94827204391f3..0000000000000 --- a/DataFormats/HGCalRecHit/src/classes.h +++ /dev/null @@ -1,3 +0,0 @@ -#include "DataFormats/Common/interface/Wrapper.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitHost.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitSoA.h" diff --git a/DataFormats/HGCalRecHit/src/classes_def.xml b/DataFormats/HGCalRecHit/src/classes_def.xml deleted file mode 100644 index 810d24bfc9a67..0000000000000 --- a/DataFormats/HGCalRecHit/src/classes_def.xml +++ /dev/null @@ -1,6 +0,0 @@ - - - - - - diff --git a/DataFormats/HGCalReco/interface/HGCalSoARecHits.h b/DataFormats/HGCalReco/interface/HGCalSoARecHits.h index ccd28f2b09da5..6e6d7b3d19c05 100644 --- a/DataFormats/HGCalReco/interface/HGCalSoARecHits.h +++ b/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 diff --git a/Geometry/HGCalMapping/plugins/alpaka/HGCalDenseIndexInfoESProducer.cc b/Geometry/HGCalMapping/plugins/alpaka/HGCalDenseIndexInfoESProducer.cc index 2b13ac2580279..9a48ed73cbec1 100644 --- a/Geometry/HGCalMapping/plugins/alpaka/HGCalDenseIndexInfoESProducer.cc +++ b/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 diff --git a/Geometry/HGCalMapping/plugins/alpaka/HGCalMappingCellESProducer.cc b/Geometry/HGCalMapping/plugins/alpaka/HGCalMappingCellESProducer.cc index 73badd50ef09d..771075aeb00b1 100644 --- a/Geometry/HGCalMapping/plugins/alpaka/HGCalMappingCellESProducer.cc +++ b/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); diff --git a/Geometry/HGCalMapping/test/testMappingModuleIndexer_cfg.py b/Geometry/HGCalMapping/test/testMappingModuleIndexer_cfg.py index e9ae812177a22..b7e1949b6e18d 100644 --- a/Geometry/HGCalMapping/test/testMappingModuleIndexer_cfg.py +++ b/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) diff --git a/RecoLocalCalo/Configuration/python/hgcalTestBeamLocalReco_cff.py b/RecoLocalCalo/Configuration/python/hgcalTestBeamLocalReco_cff.py new file mode 100644 index 0000000000000..5b79f5285ba88 --- /dev/null +++ b/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 diff --git a/RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h b/RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h index 94c119e596d63..89ef2fc1411ba 100644 --- a/RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h +++ b/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 #include using json = nlohmann::ordered_json; // ordered_json preserves key insertion order @@ -13,6 +14,57 @@ namespace hgcal { const std::string& firstkey, const std::vector& keys, const std::string& fname); + bool check_keys(const json& data, const std::vector& keys, const std::string& fname); + + // @short fill full SoA column with single value + template + 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 + void fill_SoA_column( + T* column_SoA, const std::vector& 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 + void fill_SoA_eigen_row(P& soa, const std::vector>& 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 diff --git a/RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h b/RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h index 83145c6651f1d..235f0baf155a1 100644 --- a/RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h +++ b/RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h @@ -12,8 +12,8 @@ // Host & devide HGCal RecHit data formats #include "DataFormats/HGCalDigi/interface/HGCalDigiHost.h" #include "DataFormats/HGCalDigi/interface/alpaka/HGCalDigiDevice.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitHost.h" -#include "DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h" +#include "DataFormats/HGCalReco/interface/HGCalSoARecHitsHostCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoARecHitsDeviceCollection.h" #include "CondFormats/HGCalObjects/interface/HGCalCalibrationParameterHost.h" #include "CondFormats/HGCalObjects/interface/alpaka/HGCalCalibrationParameterDevice.h" #include "CondFormats/HGCalObjects/interface/HGCalMappingParameterHost.h" @@ -29,18 +29,25 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { public: HGCalRecHitCalibrationAlgorithms(int n_blocks, int n_threads) : n_blocks_(n_blocks), n_threads_(n_threads) {} - HGCalRecHitDevice calibrate(Queue& queue, - HGCalDigiHost const& host_digis, - HGCalCalibParamDevice const& device_calib, - HGCalMappingCellParamDevice const& device_mapping, - HGCalDenseIndexInfoDevice const& device_index) const; + HGCalSoARecHitsDeviceCollection calibrate(Queue& queue, + int32_t* __restrict__ nsel, + int32_t* __restrict__ sidx, + HGCalDigiHost const& host_digis, + HGCalCalibParamDevice const& device_calib, + HGCalMappingModuleParamDevice const& device_mapmod, + HGCalMappingCellParamDevice const& device_mapping, + HGCalDenseIndexInfoDevice const& device_index) const; + + HGCalSoARecHitsDeviceCollection select(Queue& queue, + int const ndigis, + int32_t const* __restrict__ nsel, + int32_t const* __restrict__ sidx, + HGCalSoARecHitsDeviceCollection const& device_recHits) const; private: void print(HGCalDigiHost const& digis, int max = -1) const; void print_digi_device(HGCalDigiDevice const& digis, int max = -1) const; - void print_recHit_device(Queue& queue, - PortableHostCollection >::View const& recHits, - int max = -1) const; + void print_recHit_device(Queue& queue, HGCalSoARecHitsHostCollection::View const& recHits, int max = -1) const; int n_blocks_; int n_threads_; diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/BuildFile.xml b/RecoLocalCalo/HGCalRecAlgos/plugins/BuildFile.xml index 843fa22caf1f8..96845bd44da2f 100644 --- a/RecoLocalCalo/HGCalRecAlgos/plugins/BuildFile.xml +++ b/RecoLocalCalo/HGCalRecAlgos/plugins/BuildFile.xml @@ -17,7 +17,6 @@ - diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc index 3eb10ee525ccd..106e8cde70c00 100644 --- a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc +++ b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc @@ -16,33 +16,30 @@ 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, - HGCalRecHitDevice::View recHits, + HGCalSoARecHitsDeviceCollection::View recHits, + HGCalDigiDevice::ConstView digis, HGCalCalibParamDevice::ConstView calibs) const { 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(); - //recHits[idx].flags() = digiflags; + bool calibvalid = calibs[idx].valid(); + auto digiflags = digis[idx].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)); - recHits[idx].flags() = (!isAvailable) * hgcalrechit::HGCalRecHitFlags::EnergyInvalid + - (!isToAavailable) * hgcalrechit::HGCalRecHitFlags::TimeInvalid; + recHits[idx].flags() = + (!isAvailable) * HGCalRecHitFlags::kEnergyInvalid + (!isToAavailable) * HGCalRecHitFlags::kTimeInvalid; } } }; - // - 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, + HGCalSoARecHitsDeviceCollection::View recHits, + HGCalDigiDevice::ConstView digis, HGCalCalibParamDevice::ConstView calibs) const { auto adc_denoise = [&](uint32_t adc, uint32_t cm, uint32_t adcm1, float adc_ped, float cm_slope, float cm_ped, float bxm1_slope) { @@ -66,32 +63,33 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (digiflags != ::hgcal::DIGI_FLAG::NotAvailable) && calibvalid); bool useTOT((digi.tctp() == 3) && isAvailable); bool useADC(!useTOT && isAvailable); - recHits[idx].energy() = useADC * adc_denoise(digi.adc(), - digi.cm(), - digi.adcm1(), - calib.ADC_ped(), - calib.CM_slope(), - calib.CM_ped(), - calib.BXm1_slope()) + - useTOT * tot_linearization(digi.tot(), - calib.TOT_lin(), - calib.TOTtoADC(), - calib.TOT_ped(), - calib.TOT_P0(), - calib.TOT_P1(), - calib.TOT_P2()); - - //after denoising/linearization apply the MIP scale - recHits[idx].energy() *= calib.MIPS_scale(); + auto charge = useADC * adc_denoise(digi.adc(), + digi.cm(), + digi.adcm1(), + calib.ADC_ped(), + calib.CM_slope(), + calib.CM_ped(), + calib.BXm1_slope()) + + useTOT * tot_linearization(digi.tot(), + calib.TOT_lin(), + calib.TOTtoADC(), + calib.TOT_ped(), + calib.TOT_P0(), + calib.TOT_P1(), + calib.TOT_P2()); + + // after denoising/linearization apply the MIP & EM scale to convert to energy (GeV) + recHits[idx].mipEnergy() = charge * calib.MIPS_scale(); + recHits[idx].energy() = recHits[idx].mipEnergy() * 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, - HGCalRecHitDevice::View recHits, + HGCalSoARecHitsDeviceCollection::View recHits, + HGCalDigiDevice::ConstView digis, HGCalCalibParamDevice::ConstView calibs) const { auto toa_inl_corr = [&](uint32_t toa, hgcalrechit::Vector32f ctdc_p, hgcalrechit::Vector8f ftdc_p) { auto gray = toa / 256; @@ -117,20 +115,23 @@ 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; + recHits[idx].timeError() = 0.; } } }; + // @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, - HGCalRecHitDevice::View recHits, + HGCalSoARecHitsDeviceCollection::View recHits, + HGCalDigiDevice::ConstView digis, HGCalCalibParamDevice::ConstView calibs, HGCalMappingCellParamDevice::ConstView maps, HGCalDenseIndexInfoDevice::ConstView index) const { @@ -157,15 +158,15 @@ 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; } - recHits[idx + offset].flags() = hgcalrechit::HGCalRecHitFlags::Normal; + recHits[idx + offset].flags() = HGCalRecHitFlags::kNormal; recHits[idx + offset].time() = isToAavailable * time_average(recHits[idx + offset].time(), recHits[idx].time(), @@ -180,34 +181,106 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } }; + // add metadata to rechits + struct HGCalRecHitCalibrationKernel_metaData { + ALPAKA_FN_ACC void operator()(Acc1D const& acc, + HGCalSoARecHitsDeviceCollection::View recHits, + HGCalCalibParamDevice::ConstView calibs, + HGCalMappingModuleParamDevice::ConstView mapmod, + HGCalDenseIndexInfoDevice::ConstView index) const { + for (auto idx : uniform_elements(acc, recHits.metadata().size())) { + auto denseIdx = index[idx]; + auto calib = calibs[idx]; + auto recHit = recHits[idx]; + if (mapmod[denseIdx.modInfoIdx()].isSiPM()) { + recHit.dim1() = denseIdx.eta(); + recHit.dim2() = denseIdx.phi(); + } // else, isSilicon == true and eta phi values will not be used + else { + recHit.dim1() = denseIdx.x(); + recHit.dim2() = denseIdx.y(); + } + recHit.dim3() = denseIdx.z(); + recHit.sigmaNoise() = calib.Noise(); + recHit.recHitIndex() = idx; + recHit.detid() = denseIdx.detid(); + recHit.layer() = denseIdx.layer(); + } + } + }; + + // @short print RecHits for debugging struct HGCalRecHitCalibrationKernel_printRecHits { - ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalRecHitDevice::ConstView view, int size) const { + ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalSoARecHitsDeviceCollection::ConstView view, int size) const { for (int i = 0; i < size; ++i) { auto const& recHit = view[i]; - printf("%d\t%f\t%f\t%d\n", i, recHit.energy(), recHit.time(), recHit.flags()); + printf( + "i=%d\t dim1=%f\t dim2=%f\t dim3=%f\t sigmaNoise=%f\t layer=%d\t detid=%u\t mipEnergy=%f\t " + "energy=%f\t time=%f\t timeError=%f\t flags=%d\n", + i, + recHit.dim1(), + recHit.dim2(), + recHit.dim3(), + recHit.sigmaNoise(), + recHit.layer(), + recHit.detid(), + recHit.mipEnergy(), + recHit.energy(), + recHit.time(), + recHit.timeError(), + recHit.flags()); } } }; - HGCalRecHitDevice HGCalRecHitCalibrationAlgorithms::calibrate(Queue& queue, - HGCalDigiHost const& host_digis, - HGCalCalibParamDevice const& device_calib, - HGCalMappingCellParamDevice const& device_mapping, - HGCalDenseIndexInfoDevice const& device_index) const { + // + struct HGCalRecHitCalibrationKernel_countRecHits { + ALPAKA_FN_ACC void operator()(Acc1D const& acc, + int32_t* __restrict__ nsel, + int32_t* __restrict__ sidx, + HGCalSoARecHitsDeviceCollection::ConstView recHits) const { + for (auto idx : uniform_elements(acc, recHits.metadata().size())) + if (!recHits[idx].flags() && recHits[idx].energy() > 5.) + sidx[alpaka::atomicAdd(acc, nsel, 1)] = idx; + } + }; + + // + struct HGCalRecHitCalibrationKernel_copyRecHits { + ALPAKA_FN_ACC void operator()(Acc1D const& acc, + HGCalSoARecHitsDeviceCollection::View output, + HGCalSoARecHitsDeviceCollection::ConstView recHits, + int32_t const* __restrict__ sidx) const { + for (auto idx : uniform_elements(acc, output.metadata().size())) + output[idx] = recHits[sidx[idx]]; + } + }; + + // @short apply all calibration kernels + HGCalSoARecHitsDeviceCollection HGCalRecHitCalibrationAlgorithms::calibrate( + Queue& queue, + int32_t* __restrict__ nsel, + int32_t* __restrict__ sidx, + HGCalDigiHost const& host_digis, + HGCalCalibParamDevice const& device_calib, + HGCalMappingModuleParamDevice const& device_mapmod, + 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; - HGCalDigiDevice device_digis(host_digis.view().metadata().size(), queue); + auto const ndigis = host_digis.view().metadata().size(); + HGCalDigiDevice device_digis(ndigis, queue); alpaka::memcpy(queue, device_digis.buffer(), host_digis.const_buffer()); LogDebug("HGCalRecHitCalibrationAlgorithms") << "\n\nINFO -- Allocating rechits buffer and initiating values" << std::endl; - HGCalRecHitDevice device_recHits(device_digis.view().metadata().size(), queue); + HGCalSoARecHitsDeviceCollection device_recHits(ndigis, queue); // number of items per group uint32_t items = n_threads_; // use as many groups as needed to cover the whole problem - uint32_t groups = divide_up_by(device_digis.view().metadata().size(), items); + uint32_t groups = divide_up_by(ndigis, items); // map items to // - threads with a single element per thread on a GPU backend // - elements within a single thread on a CPU backend @@ -217,37 +290,75 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { alpaka::exec(queue, grid, HGCalRecHitCalibrationKernel_flagRecHits{}, - device_digis.view(), device_recHits.view(), - device_calib.view()); + device_digis.const_view(), + device_calib.const_view()); alpaka::exec(queue, grid, - HGCalRecHitCalibrationKernel_adcToCharge{}, - device_digis.view(), + HGCalRecHitCalibrationKernel_adcToEnergy{}, device_recHits.view(), - device_calib.view()); + device_digis.const_view(), + device_calib.const_view()); alpaka::exec(queue, grid, HGCalRecHitCalibrationKernel_toaToTime{}, - device_digis.view(), device_recHits.view(), - device_calib.view()); + device_digis.const_view(), + device_calib.const_view()); alpaka::exec(queue, grid, HGCalRecHitCalibrationKernel_handleCalibCell{}, - device_digis.view(), device_recHits.view(), - device_calib.view(), - device_mapping.view(), - device_index.view()); + device_digis.const_view(), + device_calib.const_view(), + device_mapping.const_view(), + device_index.const_view()); + alpaka::exec(queue, + grid, + HGCalRecHitCalibrationKernel_metaData{}, + device_recHits.view(), + device_calib.const_view(), + device_mapmod.const_view(), + device_index.const_view()); + + // select rec hits + alpaka::exec( + queue, grid, HGCalRecHitCalibrationKernel_countRecHits{}, nsel, sidx, device_recHits.const_view()); + + return device_recHits; + } + + // @short select rechits + HGCalSoARecHitsDeviceCollection HGCalRecHitCalibrationAlgorithms::select( + Queue& queue, + int const ndigis, + int32_t const* __restrict__ nsel, + int32_t const* __restrict__ sidx, + HGCalSoARecHitsDeviceCollection const& device_recHits) const { + // number of items per group + uint32_t items = n_threads_; + // use as many groups as needed to cover the whole problem + uint32_t groups = divide_up_by(ndigis, items); + // map items to + // - threads with a single element per thread on a GPU backend + // - elements within a single thread on a CPU backend + auto grid = make_workdiv(groups, items); + + HGCalSoARecHitsDeviceCollection device_selRecHits(*nsel, queue); + alpaka::exec(queue, + grid, + HGCalRecHitCalibrationKernel_copyRecHits{}, + device_selRecHits.view(), + device_recHits.const_view(), + sidx); LogDebug("HGCalRecHitCalibrationAlgorithms") << "Input recHits: " << std::endl; #ifdef EDM_ML_DEBUG int n_hits_to_print = 10; - print_recHit_device(queue, *device_recHits, n_hits_to_print); + print_recHit_device(queue, *device_selRecHits, n_hits_to_print); #endif - return device_recHits; + return device_selRecHits; } void HGCalRecHitCalibrationAlgorithms::print(HGCalDigiHost const& digis, int max) const { @@ -268,8 +379,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } } - void HGCalRecHitCalibrationAlgorithms::print_recHit_device( - Queue& queue, PortableHostCollection >::View const& recHits, int max) const { + void HGCalRecHitCalibrationAlgorithms::print_recHit_device(Queue& queue, + HGCalSoARecHitsHostCollection::View const& recHits, + int max) const { auto grid = make_workdiv(1, 1); auto size = max > 0 ? max : recHits.metadata().size(); alpaka::exec(queue, grid, HGCalRecHitCalibrationKernel_printRecHits{}, recHits, size); diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc index 821418e2d4950..597022f1ab379 100644 --- a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc +++ b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc @@ -1,7 +1,11 @@ +// Author: Izaak Neutelings (March 2024) + +// includes for CMSSW #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/ParameterSet/interface/FileInPath.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" +// includes for Alpaka #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h" #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESProducer.h" #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ModuleFactory.h" @@ -9,97 +13,120 @@ #include "HeterogeneousCore/AlpakaInterface/interface/host.h" #include "HeterogeneousCore/AlpakaInterface/interface/memory.h" -#include "DataFormats/HGCalDigi/interface/HGCalElectronicsId.h" +// includes for HGCal, calibration, and configuration parameters #include "CondFormats/HGCalObjects/interface/HGCalMappingModuleIndexer.h" #include "CondFormats/HGCalObjects/interface/HGCalCalibrationParameterHost.h" +#include "CondFormats/HGCalObjects/interface/HGCalMappingParameterHost.h" #include "CondFormats/HGCalObjects/interface/alpaka/HGCalCalibrationParameterDevice.h" #include "CondFormats/DataRecord/interface/HGCalElectronicsMappingRcd.h" #include "CondFormats/DataRecord/interface/HGCalModuleConfigurationRcd.h" // depends on HGCalElectronicsMappingRcd +#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h" // for HGCSiliconDetId::waferType #include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h" // for json, search_modkey +// includes for standard libraries #include -//#include -//#include -#include // needed to read json file with std::ifstream +#include // needed to read json file with std::ifstream +#include // for std::fill namespace ALPAKA_ACCELERATOR_NAMESPACE { namespace hgcalrechit { + using namespace ::hgcal; // for check_keys, fill_SoA class HGCalCalibrationESProducer : public ESProducer { public: HGCalCalibrationESProducer(const edm::ParameterSet& iConfig) - : ESProducer(iConfig), filename_(iConfig.getParameter("filename")) { + : ESProducer(iConfig), + filename_(iConfig.getParameter("filename")), + filenameEnergy_(iConfig.getParameter("filenameEnergyLoss")) { auto cc = setWhatProduced(this); indexToken_ = cc.consumes(iConfig.getParameter("indexSource")); + mapToken_ = cc.consumes(iConfig.getParameter("mapSource")); } static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add("filename")->setComment("Path to JSON file with calibration parameters"); + desc.add("filenameEnergyLoss") + ->setComment("Path to JSON file with energy loss & thickness corrections"); desc.add("indexSource", edm::ESInputTag("")) ->setComment("Label for module indexer to set SoA size"); + desc.add("mapSource", edm::ESInputTag("")) + ->setComment("Label for SoA with module mapper/information"); descriptions.addWithDefaultLabel(desc); } - // @short fill SoA column with data from vector for any type with some offset - template - static void fill_SoA_column( - T* column_SoA, const std::vector& 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) { - throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer") - << " Expected " << nrows << " rows, but got " << nrows_vals << "!"; - } - } else if (nrows_vals < arr_offset + nrows) { - throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer") - << " Tried to copy " << nrows << " rows to SoA with offset " << arr_offset << ", but only have " - << nrows_vals << " values in JSON!"; + // @short compute thickness correction to energy loss + // 0: CE_E_120um, 1: CE_E_200um, 2: CE_E_300um, + // 3: CE_H_120um, 4: CE_H_200um, 5: CE_H_300um + float getThicknessCorrection(const std::vector& sfs, + const uint32_t& idetid, + const int& celltype, + const std::string& fname) { + using waferType = HGCSiliconDetId::waferType; + const HGCSiliconDetId detid(idetid); + const bool isCEE = (detid.det() == DetId::HGCalEE); //layer<=17; + uint32_t idx = -1; + if (celltype == waferType::HGCalHD120) + idx = (isCEE ? 0 : 3); + else if (celltype == waferType::HGCalHD200 or celltype == waferType::HGCalLD200) + idx = (isCEE ? 1 : 4); + else if (celltype == waferType::HGCalLD300) + idx = (isCEE ? 2 : 5); + else { + cms::Exception ex("InvalidData"); + ex << "Could not find thickness correction for celltype " << celltype << " in layer" << detid.layer() + << "in '" << fname << "'!"; + ex.addContext("Calling hgcal::getThicknessCorrection()"); } - 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 - static void fill_SoA_eigen_row(P& soa, const std::vector>& values, const size_t row) { - if (row >= values.size()) - throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer") - << " Tried to copy row " << row << " to SoA, but only have " << values.size() << " values in JSON!"; - if (!values.empty() && int(values[row].size()) != soa.size()) - throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer") - << " Expected " << soa.size() << " elements in Eigen vector, but got " << values[row].size() << "!"; - for (int i = 0; i < soa.size(); i++) - soa(i) = values[row][i]; + if (idx >= sfs.size()) { + cms::Exception ex("InvalidData"); + ex << "The index of the thickness correction ()" << idx << ") for celltype " << celltype << " in layer" + << detid.layer() << "is too large for '" << fname << "'(" << sfs.size() << ")!"; + ex.addContext("Calling hgcal::getThicknessCorrection()"); + } + return sfs[idx]; } // @short create the ESProducer product: a SoA with channel-level calibration constants std::optional produce(const HGCalModuleConfigurationRcd& iRecord) { - auto const& moduleMap = iRecord.get(indexToken_); + auto const& moduleIndexer = iRecord.get(indexToken_); + auto const& moduleMapper = iRecord.get(mapToken_); edm::LogInfo("HGCalCalibrationESProducer") << "produce: filename=" << filename_.fullPath().c_str(); // load dense indexing - const uint32_t nchans = moduleMap.getMaxDataSize(); // channel-level size + const uint32_t nchans = moduleIndexer.getMaxDataSize(); // channel-level size hgcalrechit::HGCalCalibParamHost product(nchans, cms::alpakatools::host()); // load calib parameters from JSON std::ifstream infile(filename_.fullPath().c_str()); + std::ifstream infileEnergy(filenameEnergy_.fullPath().c_str()); json calib_data = json::parse(infile, nullptr, true, /*ignore_comments*/ true); - for (const auto& it : moduleMap.getTypecodeMap()) { // loop over all module typecodes - std::string const& module = it.first; // module typecode, e.g. "ML-F3PT-TX-0003" + json energy_data = json::parse(infileEnergy, nullptr, true, /*ignore_comments*/ true); + + // check keys + const std::vector energy_keys = {"dEdx", "SF_thickness_Si", "SF_thickness_SiPM"}; + check_keys(energy_data, energy_keys, filenameEnergy_.fullPath()); + const float nlayers = energy_data["dEdx"].size(); // number of absorber layers + if (nlayers != 47) // TODO: retrieve from nlayers from Geometry + edm::LogError("HGCalCalibrationESProducer") + << "Expected 47 layers, but got " << nlayers << " in " << filenameEnergy_.fullPath(); + const std::vector energylosses = energy_data["dEdx"].get>(); + + // loop over all module typecodes, e.g. "ML-F3PT-TX-0003" + for (const auto& [module, ids] : moduleIndexer.getTypecodeMap()) { + const auto [fedid, modid] = ids; // retrieve matching key (glob patterns allowed) - const auto modkey = hgcal::search_modkey(module, calib_data, filename_.fullPath()); + const auto modkey = search_modkey(module, calib_data, filename_.fullPath()); auto calib_data_ = calib_data[modkey]; // get dimensions const auto firstkey = calib_data_.begin().key(); - const uint32_t offset = moduleMap.getIndexForModuleData(module); // first channel index - const uint32_t nchans = moduleMap.getNumChannels(module); // number of channels in mapper - uint32_t nrows = calib_data_[firstkey].size(); // number of channels in JSON + const uint32_t imod = moduleIndexer.getIndexForModule(fedid, modid); // dense index in module SoA + const uint32_t offset = moduleIndexer.getIndexForModuleData(module); // first channel index + const uint32_t nchans = moduleIndexer.getNumChannels(module); // number of channels in mapper + uint32_t nrows = calib_data_[firstkey].size(); // number of channels in JSON // check number of channels & ROCs make sense if (nrows % 37 != 0) { @@ -144,14 +171,38 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { fill_SoA_eigen_row(vi.TOA_FTDC(), calib_data_["TOA_FTDC"], n); fill_SoA_eigen_row(vi.TOA_TW(), calib_data_["TOA_TW"], n); } - } + + // average energy loss of absorption layers that sandwich the sensors (do not average last layer) + // https://twiki.cern.ch/twiki/pub/CMS/HGCALSimulationAndPerformance/CalibratedRecHits.pdf + const int layer = moduleMapper.view().plane()[imod]; // counts from 1 + float dEdx = + (layer < nlayers ? (energylosses[layer - 1] + energylosses[layer]) / 2 : energylosses[nlayers - 1]); + + // compute thickness correction + float sf; + const bool isSiPM = moduleMapper.view().isSiPM()[imod]; + const int celltype = moduleMapper.view().celltype()[imod]; + const uint32_t detid = moduleMapper.view().detid()[imod]; + if (isSiPM) // scintillator + sf = energy_data["SF_thickness_SiPM"][0]; + else // Si module + sf = getThicknessCorrection(energy_data["SF_thickness_Si"], detid, celltype, filenameEnergy_.fullPath()); + edm::LogInfo("HGCalCalibrationESProducer") + << "layer=" << layer << ", celltype=" << celltype << ", isSiPM=" << isSiPM << ", dEdx=" << dEdx + << ", sf=" << sf << std::endl; + dEdx *= sf * 1e3; // apply correction and convert from MeV to GeV + fill_SoA_column_single(product.view().EM_scale(), dEdx, offset, nrows); + + } // end of loop over modules return product; } // end of produce() private: edm::ESGetToken indexToken_; + edm::ESGetToken mapToken_; const edm::FileInPath filename_; + const edm::FileInPath filenameEnergy_; }; } // namespace hgcalrechit diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitProducers.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitProducers.cc index ec91965f8381f..ad1321069afb7 100644 --- a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitProducers.cc +++ b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitProducers.cc @@ -12,7 +12,7 @@ #include // Alpaka imports -#include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/EDProducer.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/SynchronizingEDProducer.h" #include "HeterogeneousCore/AlpakaCore/interface/alpaka/EDPutToken.h" #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h" #include "HeterogeneousCore/AlpakaCore/interface/alpaka/Event.h" @@ -22,14 +22,16 @@ // includes for data formats #include "DataFormats/HGCalDigi/interface/HGCalDigiHost.h" #include "DataFormats/HGCalDigi/interface/alpaka/HGCalDigiDevice.h" -#include "DataFormats/HGCalRecHit/interface/HGCalRecHitHost.h" -#include "DataFormats/HGCalRecHit/interface/alpaka/HGCalRecHitDevice.h" +#include "DataFormats/HGCalReco/interface/HGCalSoARecHitsHostCollection.h" +#include "DataFormats/HGCalReco/interface/alpaka/HGCalSoARecHitsDeviceCollection.h" // includes for size, calibration, and configuration parameters #include "CondFormats/DataRecord/interface/HGCalElectronicsMappingRcd.h" #include "CondFormats/DataRecord/interface/HGCalModuleConfigurationRcd.h" +#include "CondFormats/DataRecord/interface/HGCalElectronicsMappingRcd.h" #include "CondFormats/HGCalObjects/interface/HGCalMappingModuleIndexer.h" #include "CondFormats/HGCalObjects/interface/HGCalCalibrationParameterHost.h" +#include "CondFormats/HGCalObjects/interface/HGCalMappingParameterHost.h" #include "CondFormats/HGCalObjects/interface/alpaka/HGCalCalibrationParameterDevice.h" #include "RecoLocalCalo/HGCalRecAlgos/interface/alpaka/HGCalRecHitCalibrationAlgorithms.h" #include "CondFormats/HGCalObjects/interface/alpaka/HGCalMappingParameterDevice.h" @@ -42,32 +44,40 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { using namespace cms::alpakatools; - class HGCalRecHitsProducer : public stream::EDProducer<> { + class HGCalRecHitsProducer : public stream::SynchronizingEDProducer<> { public: explicit HGCalRecHitsProducer(const edm::ParameterSet&); static void fillDescriptions(edm::ConfigurationDescriptions&); private: + void acquire(device::Event const&, device::EventSetup const&) override; void produce(device::Event&, device::EventSetup const&) override; edm::ESWatcher calibWatcher_; const edm::EDGetTokenT digisToken_; - edm::ESGetToken calibToken_; - device::ESGetToken mappingToken_; - device::ESGetToken indexingToken_; - const device::EDPutToken recHitsToken_; + const edm::ESGetToken calibToken_; + const device::ESGetToken mappingToken_; + const device::ESGetToken indexingToken_; + const device::ESGetToken moduleToken_; + const device::EDPutToken recHitsToken_; const HGCalRecHitCalibrationAlgorithms calibrator_; const int n_hits_scale; + int ndigis_; + cms::alpakatools::host_buffer nsel_; + std::optional> sidx_; + std::optional recHits_; }; HGCalRecHitsProducer::HGCalRecHitsProducer(const edm::ParameterSet& iConfig) - : EDProducer(iConfig), + : SynchronizingEDProducer(iConfig), digisToken_{consumes(iConfig.getParameter("digis"))}, calibToken_{esConsumes(iConfig.getParameter("calibSource"))}, mappingToken_{esConsumes(iConfig.getParameter("mappingSource"))}, indexingToken_{esConsumes(iConfig.getParameter("indexingSource"))}, + moduleToken_{esConsumes()}, recHitsToken_{produces()}, calibrator_{iConfig.getParameter("n_blocks"), iConfig.getParameter("n_threads")}, - n_hits_scale{iConfig.getParameter("n_hits_scale")} { + n_hits_scale{iConfig.getParameter("n_hits_scale")}, + nsel_{cms::alpakatools::make_host_buffer()} { #ifndef HGCAL_PERF_TEST if (n_hits_scale > 1) { throw cms::Exception("RuntimeError") << "Build with `HGCAL_PERF_TEST` flag to activate `n_hits_scale`."; @@ -87,13 +97,14 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { descriptions.addWithDefaultLabel(desc); } - void HGCalRecHitsProducer::produce(device::Event& iEvent, device::EventSetup const& iSetup) { + void HGCalRecHitsProducer::acquire(device::Event const& iEvent, device::EventSetup const& iSetup) { auto& queue = iEvent.queue(); // Read digis auto const& hostCalibParamProvider = iSetup.getData(calibToken_); auto const& deviceMappingCellParamProvider = iSetup.getData(mappingToken_); auto const& deviceIndexingParamProvider = iSetup.getData(indexingToken_); + auto const& deviceModuleInfoProvider = iSetup.getData(moduleToken_); auto const& hostDigisIn = iEvent.get(digisToken_); //printout new conditions if available @@ -159,11 +170,20 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { #ifdef HGCAL_PERF_TEST auto tmpRecHits = calibrator_.calibrate(queue, hostDigis, deviceCalibParam); - HGCalRecHitDevice recHits(oldSize, queue); + recHits_ = HGCalSoARecHitsDeviceCollection(oldSize, queue); alpaka::memcpy(queue, recHits.buffer(), tmpRecHits.const_buffer(), oldSize); #else - auto recHits = calibrator_.calibrate( - queue, hostDigis, deviceCalibParam, deviceMappingCellParamProvider, deviceIndexingParamProvider); + *nsel_ = 0; + ndigis_ = hostDigis.view().metadata().size(); + sidx_ = cms::alpakatools::make_device_buffer(queue, ndigis_); + recHits_ = calibrator_.calibrate(queue, + nsel_.data(), + sidx_->data(), + hostDigis, + deviceCalibParam, + deviceModuleInfoProvider, + deviceMappingCellParamProvider, + deviceIndexingParamProvider); #endif #ifdef EDM_ML_DEBUG @@ -173,6 +193,11 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { LogDebug("HGCalRecHitsProducer") << "Time spent calibrating " << hostDigis.view().metadata().size() << " digis: " << elapsed.count(); //<< std::endl; #endif + } + + void HGCalRecHitsProducer::produce(device::Event& iEvent, device::EventSetup const& iSetup) { + auto recHits = calibrator_.select(iEvent.queue(), ndigis_, nsel_.data(), sidx_->data(), *recHits_); + sidx_.reset(); LogDebug("HGCalRecHitsProducer") << "\n\nINFO -- storing rec hits in the event"; //<< std::endl; iEvent.emplace(recHitsToken_, std::move(recHits)); diff --git a/RecoLocalCalo/HGCalRecAlgos/src/HGCalESProducerTools.cc b/RecoLocalCalo/HGCalRecAlgos/src/HGCalESProducerTools.cc index ebfb5c106003b..54c478e8fb5e1 100644 --- a/RecoLocalCalo/HGCalRecAlgos/src/HGCalESProducerTools.cc +++ b/RecoLocalCalo/HGCalRecAlgos/src/HGCalESProducerTools.cc @@ -98,4 +98,16 @@ namespace hgcal { return iscomplete; } + // @short check if JSON data contains key + bool check_keys(const json& data, const std::vector& keys, const std::string& fname) { + bool iscomplete = true; + for (auto const& key : keys) { + if (not data.contains(key)) { + edm::LogWarning("checkkeys") << " JSON is missing key '" << key << "'! Please check file " << fname; + iscomplete = false; + } + } + return iscomplete; + } + } // namespace hgcal diff --git a/RecoLocalCalo/HGCalRecAlgos/test/BuildFile.xml b/RecoLocalCalo/HGCalRecAlgos/test/BuildFile.xml index ee8160d152499..399ef90e2342c 100644 --- a/RecoLocalCalo/HGCalRecAlgos/test/BuildFile.xml +++ b/RecoLocalCalo/HGCalRecAlgos/test/BuildFile.xml @@ -5,7 +5,6 @@ - diff --git a/RecoLocalCalo/HGCalRecAlgos/test/alpaka/HGCalRecHitESProducersTest.cc b/RecoLocalCalo/HGCalRecAlgos/test/alpaka/HGCalRecHitESProducersTest.cc index 9dd12147cd8df..656b89ef66d84 100644 --- a/RecoLocalCalo/HGCalRecAlgos/test/alpaka/HGCalRecHitESProducersTest.cc +++ b/RecoLocalCalo/HGCalRecAlgos/test/alpaka/HGCalRecHitESProducersTest.cc @@ -41,7 +41,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { static void fillDescriptions(edm::ConfigurationDescriptions&); private: - const int maxchans_, maxfeds_; + const int maxchans_, maxmods_, maxfeds_; const std::string fedjson_; // JSON file of FED configuration void produce(device::Event&, device::EventSetup const&) override; void beginRun(edm::Run const&, edm::EventSetup const&) override; @@ -54,6 +54,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { HGCalRecHitESProducersTest::HGCalRecHitESProducersTest(const edm::ParameterSet& iConfig) : EDProducer(iConfig), maxchans_(iConfig.getParameter("maxchans")), + maxmods_(iConfig.getParameter("maxmods")), maxfeds_(iConfig.getParameter("maxfeds")), fedjson_(iConfig.getParameter("fedjson")) { std::cout << "HGCalRecHitESProducersTest::HGCalRecHitESProducersTest" << std::endl; @@ -67,7 +68,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { desc.add("indexSource", edm::ESInputTag{})->setComment("Label for module indexer to set SoA size"); desc.add("configSource", edm::ESInputTag{})->setComment("Label for HGCal configuration for unpacking raw data"); desc.add("calibParamSource", edm::ESInputTag{})->setComment("Label for calibration parameters"); - desc.add("maxchans", 200)->setComment("Maximum number of channels to print"); + desc.add("maxchans", 500)->setComment("Maximum number of channels to print"); + desc.add("maxmods", 8)->setComment("Maximum number of modules to print"); desc.add("maxfeds", 25)->setComment("Maximum number of FED IDs to test"); desc.add("fedjson", "")->setComment("JSON file with FED configuration parameters"); descriptions.addWithDefaultLabel(desc); @@ -83,10 +85,18 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { return stream.str(); } + using CalibParamView = hgcalrechit::HGCalCalibParamSoALayout<>::ConstViewTemplateFreeParams<128, false, true, true>; + static void printCalPars(const CalibParamView& calibView, const int idx) { + const auto& calib = calibView[idx]; + std::cout << std::setw(6) << idx << std::setw(7) << int2hex(idx) << std::dec << std::setw(12) << calib.ADC_ped() + << std::setw(11) << calib.CM_slope() << std::setw(9) << calib.CM_ped() << std::setw(11) + << calib.EM_scale() << std::endl; + } + void HGCalRecHitESProducersTest::produce(device::Event& iEvent, device::EventSetup const& iSetup) { std::cout << "HGCalRecHitESProducersTest::produce" << std::endl; auto queue = iEvent.queue(); - auto const& moduleMap = iSetup.getData(indexerToken_); + auto const& moduleIndexer = iSetup.getData(indexerToken_); auto const& config = iSetup.getData(configToken_); // HGCalConfiguration auto const& calibParamDevice = iSetup.getData(calibParamToken_); //printf("HGCalRecHitESProducersTest::produce: time to load calibParamDevice from calib ESProducers: %f seconds\n", duration(start,now())); @@ -95,8 +105,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { std::string line = "HGCalRecHitESProducersTest::produce " + std::string(90, '-'); if (configWatcher_.check(iSetup)) { std::cout << line << std::endl; - std::cout << "HGCalRecHitESProducersTest::produce: moduleMap.getMaxDataSize()=" << moduleMap.getMaxDataSize() - << ", moduleMap.getMaxERxSize()=" << moduleMap.getMaxERxSize() << std::endl; + std::cout << "HGCalRecHitESProducersTest::produce: moduleIndexer.getMaxDataSize()=" + << moduleIndexer.getMaxDataSize() << ", moduleIndexer.getMaxERxSize()=" << moduleIndexer.getMaxERxSize() + << std::endl; // ESProducer for global HGCal configuration (structs) with header markers, etc. auto nfeds = config.feds.size(); // number of FEDs @@ -118,21 +129,36 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } } - // Alpaka ESProducer for SoA with calibration parameters with pedestals, etc. + // Alpaka ESProducer for SoA with calibration parameters with pedestals, etc. per channel std::cout << line << std::endl; int size = calibParamDevice.view().metadata().size(); std::cout << "HGCalRecHitESProducersTest::produce: device size=" << size << std::endl; - std::cout << "HGCalRecHitESProducersTest::produce: calibration constants:" << std::endl; - std::cout << " idx hex ADC_ped CM_slope CM_ped BXm1_slope" << std::endl; + std::cout << "HGCalRecHitESProducersTest::produce: calibration constants per channel:" << std::endl; + std::cout << " idx hex ADC_ped CM_slope CM_ped EM_scale" << std::endl; for (int idx = 0; idx < size; idx++) { if (idx >= maxchans_) break; - std::cout << std::setw(6) << idx << std::setw(7) << int2hex(idx) << std::dec << std::setw(12) - << calibParamDevice.view()[idx].ADC_ped() << std::setw(11) << calibParamDevice.view()[idx].CM_slope() - << std::setw(9) << calibParamDevice.view()[idx].CM_ped() << std::setw(13) - << calibParamDevice.view()[idx].BXm1_slope() << std::endl; + printCalPars(calibParamDevice.view(), idx); } - } + + // per module + int imod = 0; + std::cout << "HGCalRecHitESProducersTest::produce: calibration constants per module:" << std::endl; + std::cout << " imod typecode idx hex ADC_ped CM_slope CM_ped EM_scale" << std::endl; + for (const auto& [typecode, ids] : moduleIndexer.getTypecodeMap()) { + const auto [fedid, modid] = ids; + if (imod >= maxmods_) + break; + const uint32_t offset = moduleIndexer.getIndexForModuleData(fedid, modid, 0, 0); + uint32_t minoffset = (offset > 0 ? offset - 1 : offset); + for (uint32_t idx = minoffset; idx <= offset + 1; idx++) { + const std::string typecode_ = (idx == offset ? typecode : ""); + std::cout << std::setw(6) << imod << std::setw(18) << typecode_; + printCalPars(calibParamDevice.view(), idx); + } + imod += 1; + } + } // end if for configWatcher_::check // test JSON parser json fed_data; @@ -162,9 +188,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { const auto fedkey = hgcal::search_fedkey(fedid, fed_data, fedjson_); // search matching key fedkeys.push_back(fedkey); } - std::cout << " fedid fedkey" << std::endl; + std::cout << "HGCalRecHitESProducersTest::produce: results from search_fedkey:" << fedjson_ << std::endl; + std::cout << " fedid matched fedkey" << std::endl; for (int fedid = 0; fedid <= maxfeds_; fedid++) { - std::cout << std::setw(8) << fedid << " '" + fedkeys[fedid] + "'" << std::endl; + std::cout << std::setw(8) << fedid << " '" << fedkeys[fedid] << "'" << std::endl; } std::cout << line << std::endl; diff --git a/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py b/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py index 01bc325e696e7..ed6ccd6279543 100644 --- a/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py +++ b/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py @@ -4,6 +4,8 @@ # cmsrel CMSSW_15_1_0_pre1 # cd CMSSW_15_1_0_pre1/src/ # cmsenv +# #git cms-addpkg Geometry/HGCalMapping +# #git clone -b CalibSurrOffsetMap git@github.com:RSalvatico/Geometry-HGCalMapping.git Geometry/HGCalMapping/data # git clone https://gitlab.cern.ch/hgcal-dpg/hgcal-comm.git HGCalCommissioning # scram b -j8 # cmsRun $CMSSW_BASE/src/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py @@ -17,13 +19,14 @@ # USER OPTIONS from FWCore.ParameterSet.VarParsing import VarParsing -modmapdir = os.path.join(os.environ.get('CMSSW_BASE',''),"src/HGCalCommissioning/Configuration/data") configdir = "/eos/cms/store/group/dpg_hgcal/tb_hgcal/DPG/calibrations/SepTB2024" options = VarParsing('standard') options.register('geometry', 'ExtendedRun4D104', VarParsing.multiplicity.singleton, VarParsing.varType.string, info="geometry to use") -options.register('maxchans', 200, mytype=VarParsing.varType.int, +options.register('maxchans', 500, mytype=VarParsing.varType.int, info="maximum number of channels to print out") +options.register('maxmods', 8, mytype=VarParsing.varType.int, + info="maximum number of modules to print out") options.register('maxfeds', 30, mytype=VarParsing.varType.int, info="maximum number of FED IDs to test") options.register('fedconfig', f"{configdir}/config/config_feds_hackathon.json", mytype=VarParsing.varType.int, @@ -35,21 +38,27 @@ f"{configdir}/level0_calib_hackathon.json", mytype=VarParsing.varType.string, info="Path to calibration parameters (JSON format)") +options.register('energyloss', + f"{configdir}/../EnergyLoss/hgcal_energyloss_v16.json", + mytype=VarParsing.varType.string, + info="Path to calibration parameters (JSON format)") options.register('modules', # see https://github.com/cms-data/Geometry-HGCalMapping # or https://gitlab.cern.ch/hgcal-dpg/hgcal-comm/-/tree/master/Configuration/data/ModuleMaps #"Geometry/HGCalMapping/data/ModuleMaps/modulelocator_test.txt", # test beam with six modules - #f"{modmapdir}/ModuleMaps/modulelocator_test_2mods.txt", # fedId=0 - f"{modmapdir}/ModuleMaps/modulelocator_Sep2024TBv2.txt", # 3 layers (9 modules) + f"Geometry/HGCalMapping/data/ModuleMaps/modulelocator_Sep2024TBv2.txt", # 3 layers (9 modules) mytype=VarParsing.varType.string, info="Path to module mapper. Absolute, or relative to CMSSW src directory") -options.register('sicells', 'Geometry/HGCalMapping/data/CellMaps/WaferCellMapTraces.txt', mytype=VarParsing.varType.string, +options.register('sicells', "", mytype=VarParsing.varType.string, 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, +options.register('sipmcells', "", mytype=VarParsing.varType.string, info="Path to SiPM-on-tile cell mapper. Absolute, or relative to CMSSW src directory") options.parseArguments() +relpath = os.path.join(os.environ.get('CMSSW_BASE',''),"src") if options.params.startswith('/eos/'): - options.params = os.path.relpath(options.params,os.path.join(os.environ.get('CMSSW_BASE',''),"src")) + options.params = os.path.relpath(options.params,relpath) +if options.energyloss.startswith('/eos/'): + options.energyloss = os.path.relpath(options.energyloss,relpath) if len(options.files)==0: options.files=['file:/eos/cms/store/group/dpg_hgcal/comm_hgcal/psilva/hackhathon/23234.103_TTbar_14TeV+2026D94Aging3000/step2.root'] #options.files=['file:/eos/cms/store/group/dpg_hgcal/comm_hgcal/psilva/hackhathon/23234.103_TTbar_14TeV+2026D94Aging3000/step2.root'] @@ -62,6 +71,7 @@ print(f">>> FED config: {options.fedconfig!r}") print(f">>> ECON-D config: {options.modconfig!r}") print(f">>> Calib params: {options.params!r}") +print(f">>> Energy loss: {options.energyloss!r}") # PROCESS from Configuration.Eras.Era_Phase2C17I13M9_cff import Phase2C17I13M9 as Era_Phase2 @@ -88,6 +98,7 @@ process.load("FWCore.MessageService.MessageLogger_cfi") #process.MessageLogger.debugModules = ['*'] #"hgCalCalibrationESProducer", "hgCalConfigurationESProducer"] process.MessageLogger.cerr.threshold = 'INFO' +process.MessageLogger.cerr.noTimeStamps = True process.MessageLogger.HGCalConfigurationESProducer = { } # enable logger process.MessageLogger.HGCalCalibrationESProducer = { } process.MessageLogger.search_modkey = { } @@ -98,11 +109,12 @@ # GEOMETRY process.load(f"Configuration.Geometry.Geometry{options.geometry}Reco_cff") process.load(f"Configuration.Geometry.Geometry{options.geometry}_cff") -#process.load('Geometry.HGCalMapping.hgCalMappingIndexESSource_cfi') # old -process.load('Geometry.HGCalMapping.hgCalMappingESProducer_cfi') -process.hgCalMappingESProducer.si = cms.FileInPath(options.sicells) -process.hgCalMappingESProducer.sipm = cms.FileInPath(options.sipmcells) -process.hgCalMappingESProducer.modules = cms.FileInPath(options.modules) + +# MAPPING & INDEXING +# https://github.com/cms-sw/cmssw/blob/master/Geometry/HGCalMapping/python/hgcalmapping_cff.py +from Geometry.HGCalMapping.hgcalmapping_cff import customise_hgcalmapper +kwargs = { k: getattr(options,k) for k in ['modules','sicells','sipmcells'] if getattr(options,k)!='' } +customise_hgcalmapper(process, **kwargs) # GLOBAL CONFIGURATION ESProducers (for unpacker) #process.load("RecoLocalCalo.HGCalRecAlgos.HGCalConfigurationESProducer") @@ -112,9 +124,7 @@ fedjson=cms.string(options.fedconfig), # JSON with FED configuration parameters modjson=cms.string(options.modconfig), # JSON with ECON-D configuration parameters #passthroughMode=cms.int32(0), # ignore mismatch - #cbHeaderMarker=cms.int32(0x7f), # capture block #cbHeaderMarker=cms.int32(0x5f), # capture block - #slinkHeaderMarker=cms.int32(0x55), # S-link #slinkHeaderMarker=cms.int32(0x2a), # S-link #econdHeaderMarker=cms.int32(0x154), # ECON-D #charMode=cms.int32(1), @@ -128,7 +138,9 @@ process.hgcalCalibParamESProducer = cms.ESProducer( # ESProducer to load calibration parameters from JSON file, like pedestals 'hgcalrechit::HGCalCalibrationESProducer@alpaka', filename=cms.FileInPath(options.params), - indexSource=cms.ESInputTag('hgCalMappingESProducer','') + filenameEnergyLoss=cms.FileInPath(options.energyloss), + indexSource=cms.ESInputTag('hgCalMappingESProducer',''), + mapSource=cms.ESInputTag('hgCalMappingModuleESProducer','') ) # MAIN PROCESS @@ -140,12 +152,16 @@ configSource=cms.ESInputTag('hgcalConfigESProducer', ''), calibParamSource=cms.ESInputTag('hgcalCalibParamESProducer', ''), maxchans=cms.int32(options.maxchans), # maximum number of channels to print out + maxmods=cms.int32(options.maxmods), # maximum number of modules to print out maxfeds=cms.int32(options.maxfeds), # maximum number of FED IDs to test #fedjson=cms.string(options.fedconfig), # JSON with FED configuration parameters fedjson=cms.string(""), # use hardcoded JSON instead ) -process.t = cms.Task(process.testHGCalRecHitESProducers) -process.p = cms.Path(process.t) +process.p = cms.Path(process.testHGCalRecHitESProducers) + +#### PRINT available records (for debugging) +###process.dumpES = cms.EDAnalyzer("PrintEventSetupContent") +###process.dump = cms.Path(process.dumpES) # OUTPUT process.output = cms.OutputModule( diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc index c23f9c8fafdd1..bd672655308a2 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersAlgoWrapper.dev.cc @@ -41,7 +41,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { inputs.dim1(), inputs.dim2(), inputs.layer(), - inputs.weight(), + inputs.energy(), inputs.sigmaNoise(), inputs.detid(), outputs.rho(), diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc index a68c66ab5b315..e1cdff5480004 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalLayerClustersSoAAlgoWrapper.dev.cc @@ -25,7 +25,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { continue; } auto clIdx = input_clusters_soa[i].clusterIndex(); - alpaka::atomicAdd(acc, &outputs[clIdx].energy(), input_rechits_soa[i].weight()); + alpaka::atomicAdd(acc, &outputs[clIdx].energy(), input_rechits_soa[i].energy()); alpaka::atomicAdd(acc, &outputs[clIdx].cells(), 1); if (input_clusters_soa[i].isSeed() == 1) { outputs[clIdx].seed() = input_rechits_soa[i].detid(); @@ -54,12 +54,12 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { continue; } - alpaka::atomicAdd(acc, &outputs_service[cluster_index].total_weight(), input_rechits_soa[hit_index].weight()); + alpaka::atomicAdd(acc, &outputs_service[cluster_index].total_weight(), input_rechits_soa[hit_index].energy()); // Read the current seed index, and the associated energy. int clusterSeed = outputs_service[cluster_index].maxEnergyIndex(); - float clusterEnergy = (clusterSeed == kInvalidIndex) ? 0.f : input_rechits_soa[clusterSeed].weight(); + float clusterEnergy = (clusterSeed == kInvalidIndex) ? 0.f : input_rechits_soa[clusterSeed].energy(); - while (input_rechits_soa[hit_index].weight() > clusterEnergy) { + while (input_rechits_soa[hit_index].energy() > clusterEnergy) { // If output_service[cluster_index].maxEnergyIndex() did not change, // store the new value and exit the loop. Otherwise return the value // that has been updated, and decide again if the maximum needs to be @@ -71,7 +71,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } else { // Update the seed index and re-read the associated energy. clusterSeed = seed; - clusterEnergy = (clusterSeed == kInvalidIndex) ? 0.f : input_rechits_soa[clusterSeed].weight(); + clusterEnergy = (clusterSeed == kInvalidIndex) ? 0.f : input_rechits_soa[clusterSeed].energy(); } } // CAS } // uniform_elements @@ -105,7 +105,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { if (std::fmaf(d1, d1, d2 * d2) > positionDeltaRho2) { continue; } - float Wi = std::max(thresholdW0 + std::log(input_rechits_soa[hit_index].weight() / + float Wi = std::max(thresholdW0 + std::log(input_rechits_soa[hit_index].energy() / outputs_service[cluster_index].total_weight()), 0.f); alpaka::atomicAdd(acc, &outputs[cluster_index].x(), input_rechits_soa[hit_index].dim1() * Wi); diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc index c3d4a6f204cfd..598d4905be2a8 100644 --- a/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc +++ b/RecoLocalCalo/HGCalRecProducers/plugins/alpaka/HGCalSoARecHitsProducer.cc @@ -108,7 +108,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { entryInSoA.dim2() = position.y(); } entryInSoA.dim3() = position.z(); - entryInSoA.weight() = hgrh.energy(); + entryInSoA.energy() = hgrh.energy(); + entryInSoA.mipEnergy() = hgrh.energy(); // TODO: CHANGE TO MIP entryInSoA.sigmaNoise() = sigmaNoise; entryInSoA.layer() = layer; entryInSoA.recHitIndex() = i;