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