algorithms/calorimetry/CalorimeterHitsMerger.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence
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/IDDescriptor.h>
0016 #include <DD4hep/Objects.h>
0017 #include <DD4hep/Readout.h>
0018 #include <DD4hep/VolumeManager.h>
0019 #include <DDSegmentation/BitFieldCoder.h>
0020 #include <Evaluator/DD4hepUnits.h>
0021 #include <Math/GenVector/Cartesian3D.h>
0022 #include <Math/GenVector/DisplacementVector3D.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
0035 namespace eicrecon {
0036
0037 void CalorimeterHitsMerger::init() {
0038
0039     if (m_cfg.readout.empty()) {
0040         error("readoutClass is not provided, it is needed to know the fields in readout ids");
0041         return;
0042     }
0043
0044     try {
0045         auto id_desc = m_detector->readout(m_cfg.readout).idSpec();
0046         id_mask = 0;
0047         std::vector<std::pair<std::string, int>> ref_fields;
0048         for (size_t i = 0; i < m_cfg.fields.size(); ++i) {
0049             id_mask |= id_desc.field(m_cfg.fields[i])->mask();
0050             // use the provided id number to find ref cell, or use 0
0051             int ref = i < m_cfg.refs.size() ? m_cfg.refs[i] : 0;
0052             ref_fields.emplace_back(m_cfg.fields[i], ref);
0053         }
0054         ref_mask = id_desc.encode(ref_fields);
0055     } catch (...) {
0056         auto mess = fmt::format("Failed to load ID decoder for {}", m_cfg.readout);
0057         warning(mess);
0058 //        throw std::runtime_error(mess);
0059     }
0060     id_mask = ~id_mask;
0061     debug("ID mask in {:s}: {:#064b}", m_cfg.readout, id_mask);
0062 }
0063
0064 void CalorimeterHitsMerger::process(
0065       const CalorimeterHitsMerger::Input& input,
0066       const CalorimeterHitsMerger::Output& output) const {
0067
0068     const auto [in_hits] = input;
0069     auto [out_hits] = output;
0070
0071     // find the hits that belong to the same group (for merging)
0072     std::unordered_map<uint64_t, std::vector<std::size_t>> merge_map;
0073     std::size_t ix = 0;
0074     for (const auto &h : *in_hits) {
0075         uint64_t id = h.getCellID() & id_mask;
0076         merge_map[id].push_back(ix);
0077
0078         ix++;
0079     }
0080
0081     // sort hits by energy from large to small
0082     for (auto &it : merge_map) {
0083         std::sort(it.second.begin(), it.second.end(), [&](std::size_t ix1, std::size_t ix2) {
0084             return (*in_hits)[ix1].getEnergy() > (*in_hits)[ix2].getEnergy();
0085         });
0086     }
0087
0088     // reconstruct info for merged hits
0089     // dd4hep decoders
0090     auto volman = m_detector->volumeManager();
0091
0092     for (const auto &[id, ixs] : merge_map) {
0093         // reference fields id
0094         const uint64_t ref_id = id | ref_mask;
0095         // global positions
0096         const auto gpos = m_converter->position(ref_id);
0097         // local positions
0098         auto alignment = volman.lookupDetElement(ref_id).nominal();
0099         const auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
0100         debug("{}, {}", volman.lookupDetElement(ref_id).path(), volman.lookupDetector(ref_id).path());
0101         // sum energy
0102         float energy = 0.;
0103         float energyError = 0.;
0104         float time = 0;
0105         float timeError = 0;
0106         for (auto ix : ixs) {
0107             auto hit = (*in_hits)[ix];
0108             energy += hit.getEnergy();
0109             energyError += hit.getEnergyError() * hit.getEnergyError();
0110             time += hit.getTime();
0111             timeError += hit.getTimeError() * hit.getTimeError();
0112         }
0113         energyError = sqrt(energyError);
0114         time /= ixs.size();
0115         timeError = sqrt(timeError) / ixs.size();
0116
0117         const auto href = (*in_hits)[ixs.front()];
0118
0119         // create const vectors for passing to hit initializer list
0120         const decltype(edm4eic::CalorimeterHitData::position) position(
0121                 gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm
0122         );
0123         const decltype(edm4eic::CalorimeterHitData::local) local(
0124                 pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm
0125         );
0126
0127         out_hits->create(
0128                         href.getCellID(),
0129                         energy,
0130                         energyError,
0131                         time,
0132                         timeError,
0133                         position,
0134                         href.getDimension(),
0135                         href.getSector(),
0136                         href.getLayer(),
0137                         local); // Can do better here? Right now position is mapped on the central hit
0138     }
0139
0140     debug("Size before = {}, after = {}", in_hits->size(), out_hits->size());
0141 }
0142
0143 } // namespace eicrecon