From 6b1e136db82c8ea07a49b88eb612307c71cc7383 Mon Sep 17 00:00:00 2001 From: IzaakWN Date: Tue, 27 May 2025 19:17:20 +0200 Subject: [PATCH 1/3] add energy loss to CalibParams SoA; add module-level tester --- .../interface/HGCalCalibParamSoA.h | 1 + .../interface/HGCalESProducerTools.h | 52 +++++++ .../HGCalRecHitCalibrationESProducer.cc | 140 ++++++++++++------ .../HGCalRecAlgos/src/HGCalESProducerTools.cc | 12 ++ .../test/alpaka/HGCalRecHitESProducersTest.cc | 57 +++++-- .../test/testHGCalRecHitESProducers_cfg.py | 46 ++++-- 6 files changed, 234 insertions(+), 74 deletions(-) diff --git a/CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h b/CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h index c917cd4abc2df..4379234471fe1 100644 --- a/CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h +++ b/CondFormats/HGCalObjects/interface/HGCalCalibParamSoA.h @@ -33,6 +33,7 @@ namespace hgcalrechit { SOA_EIGEN_COLUMN(Vector8f, TOA_FTDC), // TOA fine TDC correction SOA_EIGEN_COLUMN(Vector3f, TOA_TW), // TOA timewalk correction SOA_COLUMN(float, MIPS_scale), // MIPS scale + SOA_COLUMN(float, EM_scale), // electromagnetic scale SOA_COLUMN(unsigned char, valid) // only 1 bit used: if false = mask dead channel ) using HGCalCalibParamSoA = HGCalCalibParamSoALayout<>; diff --git a/RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h b/RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h index 94c119e596d63..703a3dc0f1b39 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/MessageLogger/interface/MessageLogger.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/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc index 821418e2d4950..c7387fd91ce67 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,37 @@ 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; + fill_SoA_column_single(product.view().EM_scale(), dEdx * sf, 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/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/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..01d3db6470463 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 @@ -22,8 +24,10 @@ 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,6 +39,10 @@ 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 @@ -43,13 +51,16 @@ f"{modmapdir}/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 +73,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 +100,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 +111,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 +126,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 +140,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 +154,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( From f3c56d100b0df86d4dc192d0cf8a71530772e329 Mon Sep 17 00:00:00 2001 From: IzaakWN Date: Wed, 28 May 2025 18:46:40 +0200 Subject: [PATCH 2/3] apply energy conversion to RecHits --- .../HGCalRecHitCalibrationAlgorithms.dev.cc | 29 ++++++++++--------- .../HGCalRecHitCalibrationESProducer.cc | 3 +- 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc index 3eb10ee525ccd..2e57cbb0bc2c0 100644 --- a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc +++ b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc @@ -16,7 +16,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { using namespace cms::alpakatools; - // + // @short set RecHit flags based on DIGI flags (data corruption, zero-suppression) struct HGCalRecHitCalibrationKernel_flagRecHits { ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalDigiDevice::View digis, @@ -27,7 +27,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { bool calibvalid = calib.valid(); auto digi = digis[idx]; auto digiflags = digi.flags(); - //recHits[idx].flags() = digiflags; bool isAvailable((digiflags != ::hgcal::DIGI_FLAG::Invalid) && (digiflags != ::hgcal::DIGI_FLAG::NotAvailable) && calibvalid); bool isToAavailable((digiflags != ::hgcal::DIGI_FLAG::ZS_ToA) && @@ -38,8 +37,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } }; - // - struct HGCalRecHitCalibrationKernel_adcToCharge { + // @short subtract pedestals, linearize ADC & TOT to charge (fC), and convert to energy (GeV) + struct HGCalRecHitCalibrationKernel_adcToEnergy { ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalDigiDevice::View digis, HGCalRecHitDevice::View recHits, @@ -81,13 +80,13 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { calib.TOT_P1(), calib.TOT_P2()); - //after denoising/linearization apply the MIP scale - recHits[idx].energy() *= calib.MIPS_scale(); + // after denoising/linearization apply the MIP scale + recHits[idx].energy() *= calib.MIPS_scale() * calib.EM_scale(); } } }; - // + // @short apply INL and timewalk corrections, and convert TOA to time (ps) struct HGCalRecHitCalibrationKernel_toaToTime { ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalDigiDevice::View digis, @@ -117,16 +116,18 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { bool isToAavailable((digiflags != ::hgcal::DIGI_FLAG::ZS_ToA) && (digiflags != ::hgcal::DIGI_FLAG::ZS_ToA_ADCm1)); bool isGood(isAvailable && isToAavailable); - //INL correction + // INL correction auto toa = isGood * toa_inl_corr(digi.toa(), calib.TOA_CTDC(), calib.TOA_FTDC()); - //timewalk correction + // timewalk correction toa = isGood * toa_tw_corr(toa, recHits[idx].energy(), calib.TOA_TW()); - //toa to ps + // toa to ps recHits[idx].time() = toa * hgcalrechit::TOAtops; } } }; + // @short sum the energy of the RecHits in the calibration and surrounding cells + // see https://github.com/cms-sw/cmssw/pull/47829 struct HGCalRecHitCalibrationKernel_handleCalibCell { ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalDigiDevice::View digis, @@ -157,10 +158,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { auto cellIndex = index[idx].cellInfoIdx(); bool isCalibCell(maps[cellIndex].iscalib()); - int offset = maps[cellIndex].caliboffset(); //Calibration-to-surrounding cell offset + int offset = maps[cellIndex].caliboffset(); // calibration-to-surrounding cell offset bool is_surr_cell((offset != 0) && isAvailable && isCalibCell); - //Effectively operate only on the cell that surrounds the calibration cells + // effectively operate only on the cell that surrounds the calibration cells if (!is_surr_cell) { continue; } @@ -180,6 +181,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } }; + // @short print RecHits for debugging struct HGCalRecHitCalibrationKernel_printRecHits { ALPAKA_FN_ACC void operator()(Acc1D const& acc, HGCalRecHitDevice::ConstView view, int size) const { for (int i = 0; i < size; ++i) { @@ -189,6 +191,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } }; + // @short apply all calibration kernels HGCalRecHitDevice HGCalRecHitCalibrationAlgorithms::calibrate(Queue& queue, HGCalDigiHost const& host_digis, HGCalCalibParamDevice const& device_calib, @@ -222,7 +225,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { device_calib.view()); alpaka::exec(queue, grid, - HGCalRecHitCalibrationKernel_adcToCharge{}, + HGCalRecHitCalibrationKernel_adcToEnergy{}, device_digis.view(), device_recHits.view(), device_calib.view()); diff --git a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc index c7387fd91ce67..2e2ac8ffd3157 100644 --- a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc +++ b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationESProducer.cc @@ -190,7 +190,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { edm::LogInfo("HGCalCalibrationESProducer") << "layer=" << layer << ", celltype=" << celltype << ", isSiPM=" << isSiPM << ", dEdx=" << dEdx << ", sf=" << sf << std::endl; - fill_SoA_column_single(product.view().EM_scale(), dEdx * sf, offset, nrows); + 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 From 50c4dfcbd5eab31e4583fe7d3b8c8a97cc3f761b Mon Sep 17 00:00:00 2001 From: IzaakWN Date: Tue, 3 Jun 2025 17:30:22 +0200 Subject: [PATCH 3/3] tweaks --- .../test/testMappingModuleIndexer_cfg.py | 17 +++++++---------- .../HGCalRecHitCalibrationAlgorithms.dev.cc | 2 +- .../test/testHGCalRecHitESProducers_cfg.py | 2 +- 3 files changed, 9 insertions(+), 12 deletions(-) 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/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc index 2e57cbb0bc2c0..1b7fa844e3c8d 100644 --- a/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc +++ b/RecoLocalCalo/HGCalRecAlgos/plugins/alpaka/HGCalRecHitCalibrationAlgorithms.dev.cc @@ -80,7 +80,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { calib.TOT_P1(), calib.TOT_P2()); - // after denoising/linearization apply the MIP scale + // after denoising/linearization apply the MIP & EM scale to convert to energy (GeV) recHits[idx].energy() *= calib.MIPS_scale() * calib.EM_scale(); } } diff --git a/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py b/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py index 01d3db6470463..5ff1dc43fecb6 100644 --- a/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py +++ b/RecoLocalCalo/HGCalRecAlgos/test/testHGCalRecHitESProducers_cfg.py @@ -116,7 +116,7 @@ # 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) +customise_hgcalmapper(process, **kwargs) # GLOBAL CONFIGURATION ESProducers (for unpacker) #process.load("RecoLocalCalo.HGCalRecAlgos.HGCalConfigurationESProducer")