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