algorithms/calorimetry/CalorimeterHitsMerger.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2025 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence, Derek Anderson
0003
0004 /*
0005  *  An algorithm to group readout hits from a calorimeter
0006  *  Energy is summed
0007  *
0008  *  Author: Chao Peng (ANL), 03/31/2021
0009  */
0010
0011 #include "algorithms/calorimetry/CalorimeterHitsMerger.h"
0012
0013 #include <DD4hep/Alignments.h>
0014 #include <DD4hep/DetElement.h>
0015 #include <DD4hep/Objects.h>
0016 #include <DD4hep/Readout.h>
0017 #include <DD4hep/VolumeManager.h>
0018 #include <DDSegmentation/BitFieldCoder.h>
0019 #include <Evaluator/DD4hepUnits.h>
0020 #include <Math/GenVector/Cartesian3D.h>
0021 #include <Math/GenVector/DisplacementVector3D.h>
0022 #include <algorithms/service.h>
0023 #include <edm4hep/RawCalorimeterHit.h>
0024 #include <edm4hep/Vector3f.h>
0025 #include <fmt/core.h>
0026 #include <algorithm>
0027 #include <cmath>
0028 #include <cstddef>
0029 #include <gsl/pointers>
0030 #include <string>
0031 #include <unordered_map>
0032 #include <utility>
0033 #include <vector>
0034
0035 #include "algorithms/calorimetry/CalorimeterHitsMergerConfig.h"
0036 #include "services/evaluator/EvaluatorSvc.h"
0037
0038 namespace eicrecon {
0039
0040 void CalorimeterHitsMerger::init() {
0041
0042   if (m_cfg.readout.empty()) {
0043     error("readoutClass is not provided, it is needed to know the fields in readout ids");
0044     return;
0045   }
0046
0047   // split parameters into vectors of fields
0048   // and of transformations
0049   std::vector<std::string> fields;
0050   std::vector<std::string> transforms;
0051   for (const std::string& field_transform : m_cfg.fieldTransformations) {
0052
0053     const std::size_t isplit = field_transform.find_first_of(':');
0054     if (isplit == std::string::npos) {
0055       warning("transform '{}' ill-formatted. Format is <field>:<transformation>.", field_transform);
0056     }
0057
0058     fields.emplace_back(field_transform.substr(0, isplit));
0059     transforms.emplace_back(field_transform.substr(isplit + 1));
0060   }
0061
0062   // initialize descriptor + decoders
0063   // First, try and get the IDDescriptor. This will throw an exception if it fails.
0064   try {
0065     id_desc = m_detector->readout(m_cfg.readout).idSpec();
0066   } catch (...) {
0067     warning("Failed to get idSpec for {}", m_cfg.readout);
0068     return;
0069   }
0070   try {
0071     id_decoder = id_desc.decoder();
0072     for (const std::string& field : fields) {
0073       const short index [[maybe_unused]] = id_decoder->index(field);
0074     }
0075   } catch (...) {
0076     auto mess = fmt::format("Failed to load ID decoder for {}", m_cfg.readout);
0077     warning(mess);
0078     return;
0079   }
0080
0081   // lambda to translate IDDescriptor fields into function parameters
0082   std::function hit_transform = [this](const edm4eic::CalorimeterHit& hit) {
0083     std::unordered_map<std::string, double> params;
0084     for (const auto& name_field : id_desc.fields()) {
0085       params.emplace(name_field.first, name_field.second->value(hit.getCellID()));
0086       trace("{} = {}", name_field.first, name_field.second->value(hit.getCellID()));
0087     }
0088     return params;
0089   };
0090
0091   // loop through provided readout fields
0092   auto& svc          = algorithms::ServiceSvc::instance();
0093   std::size_t iField = 0;
0094   for (std::string& field : fields) {
0095
0096     // grab provided transformation and field
0097     const std::string field_transform = transforms.at(iField);
0098
0099     // set transformation for each field
0100     ref_maps[field] =
0101         svc.service<EvaluatorSvc>("EvaluatorSvc")->compile(field_transform, hit_transform);
0102     trace("{}: using transformation '{}'", field, field_transform);
0103     ++iField;
0104   } // end field loop
0105 }
0106
0107 void CalorimeterHitsMerger::process(const CalorimeterHitsMerger::Input& input,
0108                                     const CalorimeterHitsMerger::Output& output) const {
0109
0110   const auto [in_hits] = input;
0111   auto [out_hits]      = output;
0112
0113   // find the hits that belong to the same group (for merging)
0114   MergeMap merge_map;
0115   build_merge_map(in_hits, merge_map);
0116   debug("Merge map built: merging {} hits into {} merged hits", in_hits->size(), merge_map.size());
0117
0118   // sort hits by energy from large to small
0119   for (auto& it : merge_map) {
0120     std::sort(it.second.begin(), it.second.end(), [&](std::size_t ix1, std::size_t ix2) {
0121       return (*in_hits)[ix1].getEnergy() > (*in_hits)[ix2].getEnergy();
0122     });
0123   }
0124
0125   // reconstruct info for merged hits
0126   // dd4hep decoders
0127   auto volman = m_detector->volumeManager();
0128
0129   for (const auto& [id, ixs] : merge_map) {
0130     // reference fields id
0131     const uint64_t ref_id = id | ref_mask;
0132     // global positions
0133     const auto gpos = m_converter->position(ref_id);
0134     // local positions
0135     auto alignment = volman.lookupDetElement(ref_id).nominal();
0136     const auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
0137     debug("{}, {}", volman.lookupDetElement(ref_id).path(), volman.lookupDetector(ref_id).path());
0138     // sum energy
0139     float energy      = 0.;
0140     float energyError = 0.;
0141     float time        = 0;
0142     float timeError   = 0;
0143     for (auto ix : ixs) {
0144       auto hit = (*in_hits)[ix];
0145       energy += hit.getEnergy();
0146       energyError += hit.getEnergyError() * hit.getEnergyError();
0147       time += hit.getTime();
0148       timeError += hit.getTimeError() * hit.getTimeError();
0149     }
0150     energyError = sqrt(energyError);
0151     time /= ixs.size();
0152     timeError = sqrt(timeError) / ixs.size();
0153
0154     const auto href = (*in_hits)[ixs.front()];
0155
0156     // create const vectors for passing to hit initializer list
0157     const decltype(edm4eic::CalorimeterHitData::position) position(
0158         gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm);
0159     const decltype(edm4eic::CalorimeterHitData::local) local(
0160         pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm);
0161
0162     auto out_hit = out_hits->create(
0163         href.getCellID(), energy, energyError, time, timeError, position, href.getDimension(),
0164         href.getSector(), href.getLayer(),
0165         local); // Can do better here? Right now position is mapped on the central hit
0166
0167     // FIXME likely can do better, but for now
0168     // set related raw hit relation to be raw hit
0169     // of reference
0170     out_hit.setRawHit(href.getRawHit());
0171   }
0172
0173   debug("Size before = {}, after = {}", in_hits->size(), out_hits->size());
0174 }
0175
0176 void CalorimeterHitsMerger::build_merge_map(const edm4eic::CalorimeterHitCollection* in_hits,
0177                                             MergeMap& merge_map) const {
0178
0179   std::vector<RefField> ref_fields;
0180   for (std::size_t iHit = 0; const auto& hit : *in_hits) {
0181
0182     ref_fields.clear();
0183     for (const auto& name_field : id_desc.fields()) {
0184
0185       // apply mapping to field if provided,
0186       // otherwise copy value of field
0187       if (ref_maps.contains(name_field.first)) {
0188         ref_fields.emplace_back(name_field.first, ref_maps[name_field.first](hit));
0189       } else {
0190         ref_fields.emplace_back(name_field.first,
0191                                 id_decoder->get(hit.getCellID(), name_field.first));
0192       }
0193     }
0194
0195     // encode new cell ID and add hit to map
0196     const uint64_t ref_id = id_desc.encode(ref_fields);
0197     merge_map[ref_id].push_back(iHit);
0198     ++iHit;
0199
0200   } // end hit loop
0201
0202 } // end 'build_merge_map(edm4eic::CalorimeterHitsCollection*, MergeMap&)'
0203
0204 } // namespace eicrecon