src/components/CalorimeterHitsMerger.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck
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 #include <algorithm>
0011 #include <bitset>
0012 #include <tuple>
0013 #include <unordered_map>
0014
0015 #include "Gaudi/Property.h"
0016 #include "GaudiAlg/GaudiAlgorithm.h"
0017 #include "GaudiAlg/GaudiTool.h"
0018 #include "GaudiAlg/Transformer.h"
0019 #include "GaudiKernel/PhysicalConstants.h"
0020 #include "GaudiKernel/RndmGenerators.h"
0021 #include "GaudiKernel/ToolHandle.h"
0022
0023 #include "DDRec/CellIDPositionConverter.h"
0024 #include "DDRec/Surface.h"
0025 #include "DDRec/SurfaceManager.h"
0026 #include "DDSegmentation/BitFieldCoder.h"
0027
0028 #include "fmt/format.h"
0029 #include "fmt/ranges.h"
0030
0031 #include <k4FWCore/DataHandle.h>
0032 #include <k4Interface/IGeoSvc.h>
0033
0034 // Event Model related classes
0035 #include "edm4eic/CalorimeterHitCollection.h"
0036
0037 using namespace Gaudi::Units;
0038
0039 namespace Jug::Reco {
0040
0041 /** Calorimeter hits merging algorithm.
0042  *
0043  *  An algorithm to group readout hits from a calorimeter
0044  *  Energy is summed
0045  *
0046  *  \ingroup reco
0047  */
0048 class CalorimeterHitsMerger : public GaudiAlgorithm {
0049 private:
0050   Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
0051   Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
0052   // field names to generate id mask, the hits will be grouped by masking the field
0053   Gaudi::Property<std::vector<std::string>> u_fields{this, "fields", {"layer"}};
0054   // reference field numbers to locate position for each merged hits group
0055   Gaudi::Property<std::vector<int>> u_refs{this, "fieldRefNumbers", {}};
0056   DataHandle<edm4eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0057   DataHandle<edm4eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0058                                                                   this};
0059
0060   SmartIF<IGeoSvc> m_geoSvc;
0061   std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_converter;
0062
0063   uint64_t id_mask{0}, ref_mask{0};
0064
0065 public:
0066   CalorimeterHitsMerger(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0067     declareProperty("inputHitCollection", m_inputHitCollection, "");
0068     declareProperty("outputHitCollection", m_outputHitCollection, "");
0069   }
0070
0071   StatusCode initialize() override {
0072     if (GaudiAlgorithm::initialize().isFailure()) {
0073       return StatusCode::FAILURE;
0074     }
0075
0076     m_geoSvc = service(m_geoSvcName);
0077     if (!m_geoSvc) {
0078       error() << "Unable to locate Geometry Service. "
0079               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0080       return StatusCode::FAILURE;
0081     }
0082     m_converter = std::make_shared<const dd4hep::rec::CellIDPositionConverter>(*(m_geoSvc->getDetector()));
0083
0084     if (m_readout.value().empty()) {
0085       error() << "readoutClass is not provided, it is needed to know the fields in readout ids" << endmsg;
0086       return StatusCode::FAILURE;
0087     }
0088
0089     try {
0090       auto id_desc = m_geoSvc->getDetector()->readout(m_readout).idSpec();
0091       id_mask      = 0;
0092       std::vector<std::pair<std::string, int>> ref_fields;
0093       for (size_t i = 0; i < u_fields.size(); ++i) {
0094         id_mask |= id_desc.field(u_fields[i])->mask();
0095         // use the provided id number to find ref cell, or use 0
0096         int ref = i < u_refs.size() ? u_refs[i] : 0;
0097         ref_fields.emplace_back(u_fields[i], ref);
0098       }
0099       ref_mask = id_desc.encode(ref_fields);
0100       // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg;
0101     } catch (...) {
0102       error() << "Failed to load ID decoder for " << m_readout << endmsg;
0103       return StatusCode::FAILURE;
0104     }
0105     id_mask = ~id_mask;
0106     info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout.value(), id_mask) << endmsg;
0107     return StatusCode::SUCCESS;
0108   }
0109
0110   StatusCode execute() override {
0111     // input collections
0112     const auto& inputs = *m_inputHitCollection.get();
0113     // Create output collections
0114     auto& outputs = *m_outputHitCollection.createAndPut();
0115
0116     // find the hits that belong to the same group (for merging)
0117     std::unordered_map<long long, std::vector<size_t>> merge_map;
0118     std::size_t ix = 0;
0119     for (const auto& h : inputs) {
0120       int64_t id = h.getCellID() & id_mask;
0121       merge_map[id].push_back(ix);
0122
0123       ix++;
0124     }
0125
0126     // sort hits by energy from large to small
0127     for (auto &it : merge_map) {
0128         std::sort(it.second.begin(), it.second.end(), [&](std::size_t ix1, std::size_t ix2) {
0129             return inputs[ix1].getEnergy() > inputs[ix2].getEnergy();
0130         });
0131     }
0132
0133     // reconstruct info for merged hits
0134     // dd4hep decoders
0135     auto poscon = m_converter;
0136     auto volman = m_geoSvc->getDetector()->volumeManager();
0137
0138     for (auto& [id, ixs] : merge_map) {
0139       // reference fields id
0140       const uint64_t ref_id = id | ref_mask;
0141       // global positions
0142       const auto gpos = poscon->position(ref_id);
0143       // local positions
0144       auto alignment = volman.lookupDetElement(ref_id).nominal();
0145       const auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
0146       debug() << volman.lookupDetElement(ref_id).path() << ", " << volman.lookupDetector(ref_id).path() << endmsg;
0147       // sum energy
0148       float energy      = 0.;
0149       float energyError = 0.;
0150       float time        = 0;
0151       float timeError   = 0;
0152       for (auto ix : ixs) {
0153         auto hit = inputs[ix];
0154         energy += hit.getEnergy();
0155         energyError += hit.getEnergyError() * hit.getEnergyError();
0156         time += hit.getTime();
0157         timeError += hit.getTimeError() * hit.getTimeError();
0158       }
0159       energyError = sqrt(energyError);
0160       time /= ixs.size();
0161       timeError = sqrt(timeError) / ixs.size();
0162
0163       const auto& href = inputs[ixs.front()];
0164
0165       // create const vectors for passing to hit initializer list
0166       const decltype(edm4eic::CalorimeterHitData::position) position(
0167         gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm
0168       );
0169       const decltype(edm4eic::CalorimeterHitData::local) local(
0170         pos.x(), pos.y(), pos.z()
0171       );
0172
0173       outputs.create(
0174         href.getCellID(),
0175         energy,
0176         energyError,
0177         time,
0178         timeError,
0179         position,
0180         href.getDimension(),
0181         href.getSector(),
0182         href.getLayer(),
0183         local
0184       ); // Can do better here? Right now position is mapped on the central hit
0185     }
0186
0187     if (msgLevel(MSG::DEBUG)) {
0188       debug() << "Size before = " << inputs.size() << ", after = " << outputs.size() << endmsg;
0189     }
0190
0191     return StatusCode::SUCCESS;
0192   }
0193
0194 }; // class CalorimeterHitsMerger
0195
0196 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0197 DECLARE_COMPONENT(CalorimeterHitsMerger)
0198
0199 } // namespace Jug::Reco