From d96cbff18ba2b6e258791ffd6daf94d4b0562ba8 Mon Sep 17 00:00:00 2001 From: Derek M Anderson Date: Thu, 9 Jan 2025 10:54:37 -0500 Subject: [PATCH] Enable calorimeter hit merging by functions (#1668) ### Briefly, what does this PR introduce? This PR extends the functionality of the `CalorimeterHitsMerger` algorithm. Previously, reconstructed hits could only be merged across a given field of the readout of a calorimeter. This presents a challenge for calorimeters such as the Barrel HCal where 1. our readout has no segmentation beyond just eta & phi, and 2. it could be useful to study the response of the detector as a function of how many readout channels we gang together after reconstruction. This PR addresses this point by utilizing the `EvaluatorSvc` in a manner similar to the `adjacencyMatrix`, `peakNeighbourhoodMatrix` of the `CalorimeterIslandClustering` and the `sampFrac` of `CalorimeterHitReco`. Now the user has the ability to specify an (almost) arbitrarily complex transformation for a specific field of the readout via the `fieldTransformations` parameter, which defines both the field to transform and the function to map the indices of that field onto the desired reference indices. For example: ``` app->Add(new JOmniFactoryGeneratorT( "HcalBarrelMergedHits", {"HcalBarrelRecHits"}, {"HcalBarrelMergedHits"}, { .readout = "HcalBarrelHits", .fieldTransformations = {"phi:phi-(5*((phi/5)-floor(phi/5)))"} }, app // TODO: Remove me once fixed )); ``` Here, the `HcalBarrelMergedHits` collection will merge 5 hits (i.e. scintillator tiles for the BHCal) adjacent in phi into a one with the position and cellID of the 1st of the 5, and no hits will be merged along eta. The previous behavior of the algorithm can be recovered by simply specifying the index to be mapped onto. For example: ``` app->Add(new JOmniFactoryGeneratorT( "HcalEndcapNMergedHits", {"HcalEndcapNRecHits"}, {"HcalEndcapNMergedHits"}, { .readout = "HcalEndcapNHits", .fieldTransformations = {"layer:4", "slice:0"} }, app // TODO: Remove me once fixed )); ``` An example script of to change a transformation (and how to update the adjacency matrix accordingly) from the command-line is provided in the snippets repo [here](https://github.com/eic/snippets/blob/main/Calorimetery/CaloDebugTools/UtilityScripts/RunEICReconWithTileMerging.rb). ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [x] New feature (issue #1669 ) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [x] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No. ### Does this PR change default behavior? No. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Dmitry Kalinkin --- .../calorimetry/CalorimeterHitsMerger.cc | 117 ++++++++++++++---- .../calorimetry/CalorimeterHitsMerger.h | 25 +++- .../calorimetry/CalorimeterHitsMergerConfig.h | 3 +- src/detectors/BHCAL/BHCAL.cc | 18 +++ src/detectors/EHCAL/EHCAL.cc | 3 +- src/detectors/FHCAL/FHCAL.cc | 3 +- .../CalorimeterHitsMerger_factory.h | 3 +- src/services/io/podio/JEventProcessorPODIO.cc | 1 + 8 files changed, 139 insertions(+), 34 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterHitsMerger.cc b/src/algorithms/calorimetry/CalorimeterHitsMerger.cc index 0267d5f2bc..f43b75cce4 100644 --- a/src/algorithms/calorimetry/CalorimeterHitsMerger.cc +++ b/src/algorithms/calorimetry/CalorimeterHitsMerger.cc @@ -1,5 +1,5 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence +// Copyright (C) 2022 - 2025 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence, Derek Anderson /* * An algorithm to group readout hits from a calorimeter @@ -12,7 +12,6 @@ #include #include -#include #include #include #include @@ -20,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -31,6 +31,7 @@ #include #include "algorithms/calorimetry/CalorimeterHitsMergerConfig.h" +#include "services/evaluator/EvaluatorSvc.h" namespace eicrecon { @@ -41,24 +42,64 @@ void CalorimeterHitsMerger::init() { return; } + // split parameters into vectors of fields + // and of transformations + std::vector fields; + std::vector transforms; + for (const std::string& field_transform : m_cfg.fieldTransformations) { + + const std::size_t isplit = field_transform.find_first_of(':'); + if (isplit == std::string::npos) { + warning("transform '{}' ill-formatted. Format is :.", field_transform); + } + + + fields.emplace_back( + field_transform.substr(0, isplit) + ); + transforms.emplace_back( + field_transform.substr(isplit + 1) + ); + } + + + // initialize descriptor + decoders try { - auto id_desc = m_detector->readout(m_cfg.readout).idSpec(); - id_mask = 0; - std::vector> ref_fields; - for (size_t i = 0; i < m_cfg.fields.size(); ++i) { - id_mask |= id_desc.field(m_cfg.fields[i])->mask(); - // use the provided id number to find ref cell, or use 0 - int ref = i < m_cfg.refs.size() ? m_cfg.refs[i] : 0; - ref_fields.emplace_back(m_cfg.fields[i], ref); - } - ref_mask = id_desc.encode(ref_fields); + id_desc = m_detector->readout(m_cfg.readout).idSpec(); + id_decoder = id_desc.decoder(); + for (const std::string& field : fields) { + const short index = id_decoder->index(field); + } } catch (...) { - auto mess = fmt::format("Failed to load ID decoder for {}", m_cfg.readout); - warning(mess); + auto mess = fmt::format("Failed to load ID decoder for {}", m_cfg.readout); + warning(mess); // throw std::runtime_error(mess); } - id_mask = ~id_mask; - debug("ID mask in {:s}: {:#064b}", m_cfg.readout, id_mask); + + // lambda to translate IDDescriptor fields into function parameters + std::function hit_transform = [this](const edm4eic::CalorimeterHit& hit) { + std::unordered_map params; + for (const auto& name_field : id_desc.fields()) { + params.emplace(name_field.first, name_field.second->value(hit.getCellID())); + trace("{} = {}", name_field.first, name_field.second->value(hit.getCellID())); + } + return params; + }; + + // loop through provided readout fields + auto& svc = algorithms::ServiceSvc::instance(); + for (std::size_t iField = 0; std::string& field : fields) { + + // grab provided transformation and field + const std::string field_transform = transforms.at(iField); + auto name_field = id_desc.field(field); + + // set transformation for each field + ref_maps[field] = svc.service("EvaluatorSvc")->compile(field_transform, hit_transform); + trace("{}: using transformation '{}'", field, field_transform); + ++iField; + } // end field loop + } void CalorimeterHitsMerger::process( @@ -69,14 +110,9 @@ void CalorimeterHitsMerger::process( auto [out_hits] = output; // find the hits that belong to the same group (for merging) - std::unordered_map> merge_map; - std::size_t ix = 0; - for (const auto &h : *in_hits) { - uint64_t id = h.getCellID() & id_mask; - merge_map[id].push_back(ix); - - ix++; - } + MergeMap merge_map; + build_merge_map(in_hits, merge_map); + debug("Merge map built: merging {} hits into {} merged hits", in_hits->size(), merge_map.size()); // sort hits by energy from large to small for (auto &it : merge_map) { @@ -140,4 +176,37 @@ void CalorimeterHitsMerger::process( debug("Size before = {}, after = {}", in_hits->size(), out_hits->size()); } +void CalorimeterHitsMerger::build_merge_map( + const edm4eic::CalorimeterHitCollection* in_hits, + MergeMap& merge_map) const { + + std::vector ref_fields; + for (std::size_t iHit = 0; const auto& hit : *in_hits) { + + ref_fields.clear(); + for (std::size_t iField = 0; const auto& name_field : id_desc.fields()) { + + // apply mapping to field if provided, + // otherwise copy value of field + if (ref_maps.count(name_field.first) > 0) { + ref_fields.push_back( + {name_field.first, ref_maps[name_field.first](hit)} + ); + } else { + ref_fields.push_back( + {name_field.first, id_decoder->get(hit.getCellID(), name_field.first)} + ); + } + ++iField; + } + + // encode new cell ID and add hit to map + const uint64_t ref_id = id_desc.encode(ref_fields); + merge_map[ref_id].push_back(iHit); + ++iHit; + + } // end hit loop + +} // end 'build_merge_map(edm4eic::CalorimeterHitsCollection*, MergeMap&)' + } // namespace eicrecon diff --git a/src/algorithms/calorimetry/CalorimeterHitsMerger.h b/src/algorithms/calorimetry/CalorimeterHitsMerger.h index 763d9d199c..4de553e4a8 100644 --- a/src/algorithms/calorimetry/CalorimeterHitsMerger.h +++ b/src/algorithms/calorimetry/CalorimeterHitsMerger.h @@ -1,5 +1,5 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence +// Copyright (C) 2022 - 2025 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence, Derek Anderson /* * An algorithm to group readout hits from a calorimeter @@ -11,20 +11,33 @@ #pragma once #include +#include #include +#include #include #include #include #include +#include +#include #include +#include #include #include +#include +#include +#include #include "CalorimeterHitsMergerConfig.h" #include "algorithms/interfaces/WithPodConfig.h" namespace eicrecon { + // aliases for convenience + using MergeMap = std::unordered_map>; + using RefField = std::pair; + using MapFunc = std::function; + using CalorimeterHitsMergerAlgorithm = algorithms::Algorithm< algorithms::Input< edm4eic::CalorimeterHitCollection @@ -49,12 +62,20 @@ namespace eicrecon { void process(const Input&, const Output&) const final; private: - uint64_t id_mask{0}, ref_mask{0}; + uint64_t ref_mask{0}; + + private: + mutable std::map ref_maps; + dd4hep::IDDescriptor id_desc; + dd4hep::BitFieldCoder* id_decoder; private: const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()}; const dd4hep::rec::CellIDPositionConverter* m_converter{algorithms::GeoSvc::instance().cellIDPositionConverter()}; + private: + void build_merge_map(const edm4eic::CalorimeterHitCollection* in_hits, MergeMap& merge_map) const; + }; } // namespace eicrecon diff --git a/src/algorithms/calorimetry/CalorimeterHitsMergerConfig.h b/src/algorithms/calorimetry/CalorimeterHitsMergerConfig.h index 9d5f92df6c..bb644b6eac 100644 --- a/src/algorithms/calorimetry/CalorimeterHitsMergerConfig.h +++ b/src/algorithms/calorimetry/CalorimeterHitsMergerConfig.h @@ -11,8 +11,7 @@ namespace eicrecon { struct CalorimeterHitsMergerConfig { std::string readout{""}; - std::vector fields{}; - std::vector refs{}; + std::vector fieldTransformations{}; }; diff --git a/src/detectors/BHCAL/BHCAL.cc b/src/detectors/BHCAL/BHCAL.cc index 39871793fc..f17e908d8d 100644 --- a/src/detectors/BHCAL/BHCAL.cc +++ b/src/detectors/BHCAL/BHCAL.cc @@ -12,6 +12,7 @@ #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h" #include "factories/calorimetry/CalorimeterHitDigi_factory.h" #include "factories/calorimetry/CalorimeterHitReco_factory.h" +#include "factories/calorimetry/CalorimeterHitsMerger_factory.h" #include "factories/calorimetry/CalorimeterIslandCluster_factory.h" #include "factories/calorimetry/CalorimeterTruthClustering_factory.h" #include "factories/calorimetry/TrackClusterMergeSplitter_factory.h" @@ -66,6 +67,7 @@ extern "C" { }, app // TODO: Remove me once fixed )); + app->Add(new JOmniFactoryGeneratorT( "HcalBarrelRecHits", {"HcalBarrelRawHits"}, {"HcalBarrelRecHits"}, { @@ -83,10 +85,25 @@ extern "C" { }, app // TODO: Remove me once fixed )); + + // -------------------------------------------------------------------- + // If needed, merge adjacent phi tiles into towers. By default, + // NO merging will be done. This can be changed at runtime. + // -------------------------------------------------------------------- + app->Add(new JOmniFactoryGeneratorT( + "HcalBarrelMergedHits", {"HcalBarrelRecHits"}, {"HcalBarrelMergedHits"}, + { + .readout = "HcalBarrelHits", + .fieldTransformations = {"phi:phi"} + }, + app // TODO: Remove me once fixed + )); + app->Add(new JOmniFactoryGeneratorT( "HcalBarrelTruthProtoClusters", {"HcalBarrelRecHits", "HcalBarrelHits"}, {"HcalBarrelTruthProtoClusters"}, app // TODO: Remove me once fixed )); + app->Add(new JOmniFactoryGeneratorT( "HcalBarrelIslandProtoClusters", {"HcalBarrelRecHits"}, {"HcalBarrelIslandProtoClusters"}, { @@ -179,5 +196,6 @@ extern "C" { app // TODO: Remove me once fixed ) ); + } } diff --git a/src/detectors/EHCAL/EHCAL.cc b/src/detectors/EHCAL/EHCAL.cc index 49e21f6beb..7d39c01472 100644 --- a/src/detectors/EHCAL/EHCAL.cc +++ b/src/detectors/EHCAL/EHCAL.cc @@ -71,8 +71,7 @@ extern "C" { "HcalEndcapNMergedHits", {"HcalEndcapNRecHits"}, {"HcalEndcapNMergedHits"}, { .readout = "HcalEndcapNHits", - .fields = {"layer", "slice"}, - .refs = {4, 0}, // place merged hits at ~1 interaction length deep + .fieldTransformations = {"layer:4", "slice:0"} }, app // TODO: Remove me once fixed )); diff --git a/src/detectors/FHCAL/FHCAL.cc b/src/detectors/FHCAL/FHCAL.cc index 3f9f145528..36f7fb4797 100644 --- a/src/detectors/FHCAL/FHCAL.cc +++ b/src/detectors/FHCAL/FHCAL.cc @@ -83,8 +83,7 @@ extern "C" { "HcalEndcapPInsertMergedHits", {"HcalEndcapPInsertRecHits"}, {"HcalEndcapPInsertMergedHits"}, { .readout = "HcalEndcapPInsertHits", - .fields = {"layer", "slice"}, - .refs = {1, 0}, + .fieldTransformations = {"layer:1", "slice:0"} }, app // TODO: Remove me once fixed )); diff --git a/src/factories/calorimetry/CalorimeterHitsMerger_factory.h b/src/factories/calorimetry/CalorimeterHitsMerger_factory.h index 9f6aefcd1b..7105376532 100644 --- a/src/factories/calorimetry/CalorimeterHitsMerger_factory.h +++ b/src/factories/calorimetry/CalorimeterHitsMerger_factory.h @@ -20,8 +20,7 @@ class CalorimeterHitsMerger_factory : public JOmniFactory m_hits_output {this}; ParameterRef m_readout {this, "readout", config().readout}; - ParameterRef> m_fields {this, "fields", config().fields}; - ParameterRef> m_refs {this, "refs", config().refs}; + ParameterRef> m_field_transformations {this, "fieldTransformations", config().fieldTransformations}; Service m_algorithmsInit {this}; diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 2575c326a9..1d6cbf5d4e 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -290,6 +290,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "LFHCALSplitMergeClusterAssociations", "HcalBarrelRawHits", "HcalBarrelRecHits", + "HcalBarrelMergedHits", "HcalBarrelClusters", "HcalBarrelClusterAssociations", "HcalBarrelSplitMergeClusters",